/* Last edited on 2024-11-29 10:45:00 by stolfi */

void gausol_test_tools_check_extract_solution
  ( uint32_t m,
    uint32_t np,
    double AB[],
    uint32_t p,
    double X[],
    uint32_t rank_ext,
    uint32_t lead[],
    double X_ref[],
    double Bmax[]
  )
  { /* The normalization conditions should still hold: */
    check_normalize(m, np, AB, p, X_ref, Bmax);

    /* Check the {gausol_extract_solution} post-conditions: */
    double tol[p]; /* Max roundoff error for each column of {X} and {B}. */
    for (uint32_t k = 0;  k < p; k++) { tol[k] = max_residual_roundoff(m, np, Bmax[k]); }

    uint32_t n = np - p;
    uint32_t neq = 0; /* Number of non-zero equations. */
    uint32_t j = 0; /* Next candidate for pivot. */
    for (uint32_t i = 0;  i < m; i++)
      { while ((j < n) && (j < lead[i]))
          { for (uint32_t k = 0; k < p; k++)
              { double Xjk = X[j*p + k];
                if ((! isfinite(Xjk)) || (Xjk != 0.0))
                  { fprintf(stderr, "**element {X[%d,%d]} = %24.16e should be zero\n", i, j, Xij);
                    ok = FALSE;
                  }
              }
          }
        if (j < n)
          { /* {AB[i,j]} is the pivot on line {i}. */
            assert(j == lead[i]);
            double ABij = AB[i*np + jr]; /* Pivot element. */
            assert(isfinite(ABij));
            if ((! isfinite(ABij)) || (ABij == 0.0))
              { fprintf(stderr, "** pivot {A[%d,%d]}} = %24.16e is invalid or zero\n", i, j, ABij);
                ok = FALSE;
              }
            else
              { /* We have one more effective equation: */
                for (uint32_t k = 0; k < p; k++)
                  { double Aij = AB[i*np + j];
                    double Bik = AB[i*np + n + k];
                    double Xjk_ref = X_ref[j*p + k]; /* Intended solution. */
                    double Bik_exp = Aij * Xjk_ref;  /* Expected {B}. */
                    double Xjk_cmp = X[j*p + k];     /* Solution by library. */
                    double Bik_cmp = Aij * Xjk_cmp;  /* Its {B}. */
                    /* The residual must be zero or nearly so: */
                    if (fabs(Bik_exp - Bik_cmp) > tol[k])
                      { fprintf(stderr, "** A[%d,%d] = %24.16e", i, j, Aij);
                        fprintf(stderr, " B[%d,%d] = %24.16e", i, k, Bik);
                        fprintf(stderr, " (A X)[%d,%d] = %24.16e B[%d,%d] = %24.16e err = %24.16e tol = %24.16e\n", i, k, Yik, tol[k]);
                        fprintf(stderr, "** (A X)[%d,%d] = %24.16e B[%d,%d] = %24.16e err = %24.16e tol = %24.16e\n", i, k, Yik, tol[k]);
                        demand(FALSE, "** this residual should be zero");
                  }
              }
          }
        else
          { /* Equation was ignored, nothing to check: */ }

        /* In any case, the next row's pivot must be at least one col to the right: */
        if (jr < UINT32_MAX) { jr++; }
      }

    /* Check number of equations used. */
    if (neq != rank_ext)
      { fprintf(stderr, " number of useful equations = %d  returned by {extract_solution} = %d\n", neq, rank_ext);
        demand(FALSE, "** counts don't match");
      }
  }

double gausol_test_tools_max_residual_roundoff(uint32_t m, uint32_t n, double bmax)
  { /* Assumes that the magnitude of the products {A[i,j]*x[j]}
      is comparable to the magnitude {bmax} of the right-hand-side
      vector. Also assumes that roundoff errors are independent: */
    return 1.0e-14 * bmax * sqrt(n);
  }
]   
void gausol_test_tools_unpack_system
  ( uint32_t m, uint32_t n, uint32_t p,
    double M[], double A[], double B[],
    uint32_t prow[], uint32_t pcol[],
    uint32_t rank,
    char *head,
    bool_t verbose
  );
  /* Given an {m×(n+p)} matrix {M}, copies the first {m} columns into
    {A[0..m*n-1]} and the last {p} columns into {B[0..m*p-1]}. If
    {verbose} is true, also prints {A} and {B} side by side with the
    given {head} line, first as-is and then (if {prow} is not {NULL}) in
    permuted and partitioned form, according to {prow,pcol,rank}. */

void gausol_test_tools_unpack_system
  ( uint32_t m, uint32_t n, uint32_t p,
    double M[], double A[], double B[],
    uint32_t prow[], uint32_t pcol[],
    uint32_t rank,
    char *head,
    bool_t verbose
  )  
  { 
    uint32_t np = n + p;
    for (uint32_t i = 0; i < m; i++)
      { for (uint32_t j = 0; j < np; j++)
          { double *el = (j < n ? &(A[i*n + j]) : &(B[i*p + (j - n)]));
            (*el) = M[i*np + j];
          }
      }
    if (verbose) 
      { gausol_print_system
          ( stderr, 6, "%12.6f", head, 
            m,prow,rank, n,pcol,rank, "A",A,
            p, "B",B, 
            0,NULL,NULL, ""
          ); 
      }
  }

void gausol_solve_packed
  ( uint32_t m,
    uint32_t n,
    uint32_t p,
    double AB[],
    double X[],
    bool_t pivot_rows,
    bool_t pivot_cols,
    double tiny,
    uint32_t *rank_P,
    double *det_P
  );
  /* Like {gausol_solve}, except that 
  
      * The arrrays {A} and {B} are taken to be the first {n} and the
      last {p} columns of the {m×(n+p)} array {AB[0..m*(n+p)-1]}.
    
      * The contents of the array {AB} is lost (modified by the Gauss
      elimination method).
      
      * If {rank_P} is not {NULL}, returns in {*rank_P} the
      number {rank} of equations (in {0..min(m,n)})
      that were used in the computation of {X}; */
void check_extract_solution
  ( uint32_t m,
    uint32_t n,
    double M[],
    uint32_t p,
    double X[],
    uint32_t rank_ext,
    uint32_t lead[],
    double X_ref[],
    double Bmax[]
  );
  /* Checks the outcome of {gausol_extract_solution}. The parameter {X}
    should be the solution computed by {gausol_extract_solution},
    {X_ref} should be the `true' solution, and {rank} should be its
    return value. */