libgenua
Basic Geometry, Numerical Algorithms and Interfaces

Craig's method for underdetermined systems.
CRAIG finds a solution x to the linear equation Ax = b, where A is a real matrix with m rows and n columns, and b is a real mvector. If A is square and nonsingular, CRAIG finds the unique solution x = A(inv)b. If the system Ax = b is underdetermined (i.e. there are many solutions), CRAIG finds the solution of minimum Euclidean length, namely x = A inv(A A') b. Thus, CRAIG solves the problem
min x'x subject to Ax = b.
y returns a vector satisfying A'y = x. Hence AA'y = b.
A is an m by n matrix (ideally sparse), or a function handle such that y = A(x,1) returns y = A*x (where x will be an nvector); y = A(x,2) returns y = A'*x (where x will be an mvector).
CRAIG uses an iterative (conjugategradientlike) method. For further information, see
08 Apr 2003: First craig.m derived from Fortran 77 version of craig1.for. Michael Saunders, Systems Optimization Laboratory, Dept of MS&E, Stanford University. 09 Apr 2003: Experimenting on singular systems (for inverse iteration). Separately, on fullrank Ax = b, "Acond" seems to overestimate cond(A) drastically. 02 Oct 2006: Output y such that x = A'y (already in the f77 version). 15 Aug 2014: A can now be a matrix or a function handle, as in lsqrSOL. 28 Aug 2014: Fixed glitches found by Dominique Orban.
#include <craig.h>
Public Member Functions  
int  solve (const Operator &A, const DVector< Scalar > &b, DVector< Scalar > &x) 
minimize x^2 such that Ax = b  
Public Member Functions inherited from SolIterativeSolver  
void  tolerance (int maxiter, double atoler, double btoler, double conlim=1e8) 
set convergence tolerances and maximum number of iterations  
Additional Inherited Members  
Static Public Member Functions inherited from SolIterativeSolver  
static const char *  statusMessage (int code) 
access text corresponding to error code  
static bool  success (int code) 
does exit code indicate succes (in some sense, at least)?  
Protected Attributes inherited from SolIterativeSolver  
double  m_atol = 1e6 
convergence criteria  
double  m_btol = 1e6 
convergence criteria  
double  m_conlim = 1e8 
condition number limit (actually, limit for the estimated condition)  
size_t  m_maxiter = 128 
permitted number of iteration  