📄 lsqr.c
字号:
/* linalg/lsqr.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"
/* Table of constant values */
static integer c__1 = 1;
static integer c__2 = 2;
/* From arpa!sol-michael.stanford.edu!mike 5 May 89 23:53:00 PDT */
/*< >*/
/* Subroutine */ int lsqr_(integer *m, integer *n,
int (*aprod)(v3p_netlib_integer*,
v3p_netlib_integer*,
v3p_netlib_integer*,
v3p_netlib_doublereal*,
v3p_netlib_doublereal*,
v3p_netlib_integer*,
v3p_netlib_integer*,
v3p_netlib_integer*,
v3p_netlib_doublereal*,
void*),
doublereal *
damp, integer *leniw, integer *lenrw, integer *iw, doublereal *rw,
doublereal *u, doublereal *v, doublereal *w, doublereal *x,
doublereal *se, doublereal *atol, doublereal *btol, doublereal *
conlim, integer *itnlim, integer *nout, integer *istop, integer *itn,
doublereal *anorm, doublereal *acond, doublereal *rnorm, doublereal *
arnorm, doublereal *xnorm, void* userdata)
{
/* System generated locals */
integer i__1;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
integer i__;
doublereal t, z__, t1, t2, t3, cs, sn, cs1, cs2, sn1, sn2, phi, rho, tau,
psi, rhs, res1, res2, alfa, beta, zbar, ctol, rtol;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
doublereal test1, test2, test3, gamma;
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *);
doublereal delta, theta, bnorm;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *);
integer nconv, nstop;
doublereal rhbar1, rhbar2, gambar, phibar, rhobar, bbnorm, ddnorm, dampsq,
xxnorm;
/*< EXTERNAL APROD >*/
/*< INTEGER M, N, LENIW, LENRW, ITNLIM, NOUT, ISTOP, ITN >*/
/*< INTEGER IW(LENIW) >*/
/*< >*/
/* ----------------------------------------------------------------------- */
/* LSQR finds a solution x to the following problems: */
/* 1. Unsymmetric equations -- solve A*x = b */
/* 2. Linear least squares -- solve A*x = b */
/* in the least-squares sense */
/* 3. Damped least squares -- solve ( A )*x = ( b ) */
/* ( damp*I ) ( 0 ) */
/* in the least-squares sense */
/* where A is a matrix with m rows and n columns, b is an */
/* m-vector, and damp is a scalar. (All quantities are real.) */
/* The matrix A is intended to be large and sparse. It is accessed */
/* by means of subroutine calls of the form */
/* CALL APROD ( mode, m, n, x, y, LENIW, LENRW, IW, RW ) */
/* which must perform the following functions: */
/* If MODE = 1, compute y = y + A*x. */
/* If MODE = 2, compute x = x + A(transpose)*y. */
/* The vectors x and y are input parameters in both cases. */
/* If mode = 1, y should be altered without changing x. */
/* If mode = 2, x should be altered without changing y. */
/* The parameters LENIW, LENRW, IW, RW may be used for workspace */
/* as described below. */
/* The rhs vector b is input via U, and subsequently overwritten. */
/* Note: LSQR uses an iterative method to approximate the solution. */
/* The number of iterations required to reach a certain accuracy */
/* depends strongly on the scaling of the problem. Poor scaling of */
/* the rows or columns of A should therefore be avoided where */
/* possible. */
/* For example, in problem 1 the solution is unaltered by */
/* row-scaling. If a row of A is very small or large compared to */
/* the other rows of A, the corresponding row of ( A b ) should be */
/* scaled up or down. */
/* In problems 1 and 2, the solution x is easily recovered */
/* following column-scaling. Unless better information is known, */
/* the nonzero columns of A should be scaled so that they all have */
/* the same Euclidean norm (e.g., 1.0). */
/* In problem 3, there is no freedom to re-scale if damp is */
/* nonzero. However, the value of damp should be assigned only */
/* after attention has been paid to the scaling of A. */
/* The parameter damp is intended to help regularize */
/* ill-conditioned systems, by preventing the true solution from */
/* being very large. Another aid to regularization is provided by */
/* the parameter ACOND, which may be used to terminate iterations */
/* before the computed solution becomes very large. */
/* Notation */
/* -------- */
/* The following quantities are used in discussing the subroutine */
/* parameters: */
/* Abar = ( A ), bbar = ( b ) */
/* ( damp*I ) ( 0 ) */
/* r = b - A*x, rbar = bbar - Abar*x */
/* rnorm = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) */
/* = norm( rbar ) */
/* RELPR = the relative precision of floating-point arithmetic */
/* on the machine being used. For example, on the IBM 370, */
/* RELPR is about 1.0E-6 and 1.0D-16 in single and double */
/* precision respectively. */
/* LSQR minimizes the function rnorm with respect to x. */
/* Parameters */
/* ---------- */
/* M input m, the number of rows in A. */
/* N input n, the number of columns in A. */
/* APROD external See above. */
/* DAMP input The damping parameter for problem 3 above. */
/* (DAMP should be 0.0 for problems 1 and 2.) */
/* If the system A*x = b is incompatible, values */
/* of DAMP in the range 0 to sqrt(RELPR)*norm(A) */
/* will probably have a negligible effect. */
/* Larger values of DAMP will tend to decrease */
/* the norm of x and reduce the number of */
/* iterations required by LSQR. */
/* The work per iteration and the storage needed */
/* by LSQR are the same for all values of DAMP. */
/* LENIW input The length of the workspace array IW. */
/* LENRW input The length of the workspace array RW. */
/* IW workspace An integer array of length LENIW. */
/* RW workspace A real array of length LENRW. */
/* Note: LSQR does not explicitly use the previous four */
/* parameters, but passes them to subroutine APROD for */
/* possible use as workspace. If APROD does not need */
/* IW or RW, the values LENIW = 1 or LENRW = 1 should */
/* be used, and the actual parameters corresponding to */
/* IW or RW may be any convenient array of suitable type. */
/* U(M) input The rhs vector b. Beware that U is */
/* over-written by LSQR. */
/* V(N) workspace */
/* W(N) workspace */
/* X(N) output Returns the computed solution x. */
/* SE(N) output Returns standard error estimates for the */
/* components of X. For each i, SE(i) is set */
/* to the value rnorm * sqrt( sigma(i,i) / T ), */
/* where sigma(i,i) is an estimate of the i-th */
/* diagonal of the inverse of Abar(transpose)*Abar */
/* and T = 1 if m .le. n, */
/* T = m - n if m .gt. n and damp = 0, */
/* T = m if damp .ne. 0. */
/* ATOL input An estimate of the relative error in the data */
/* defining the matrix A. For example, */
/* if A is accurate to about 6 digits, set */
/* ATOL = 1.0E-6 . */
/* BTOL input An extimate of the relative error in the data */
/* defining the rhs vector b. For example, */
/* if b is accurate to about 6 digits, set */
/* BTOL = 1.0E-6 . */
/* CONLIM input An upper limit on cond(Abar), the apparent */
/* condition number of the matrix Abar. */
/* Iterations will be terminated if a computed */
/* estimate of cond(Abar) exceeds CONLIM. */
/* This is intended to prevent certain small or */
/* zero singular values of A or Abar from */
/* coming into effect and causing unwanted growth */
/* in the computed solution. */
/* CONLIM and DAMP may be used separately or */
/* together to regularize ill-conditioned systems. */
/* Normally, CONLIM should be in the range */
/* 1000 to 1/RELPR. */
/* Suggested value: */
/* CONLIM = 1/(100*RELPR) for compatible systems, */
/* CONLIM = 1/(10*sqrt(RELPR)) for least squares. */
/* Note: If the user is not concerned about the parameters */
/* ATOL, BTOL and CONLIM, any or all of them may be set */
/* to zero. The effect will be the same as the values */
/* RELPR, RELPR and 1/RELPR respectively. */
/* ITNLIM input An upper limit on the number of iterations. */
/* Suggested value: */
/* ITNLIM = n/2 for well-conditioned systems */
/* with clustered singular values, */
/* ITNLIM = 4*n otherwise. */
/* NOUT input File number for printed output. If positive, */
/* a summary will be printed on file NOUT. */
/* ISTOP output An integer giving the reason for termination: */
/* 0 x = 0 is the exact solution. */
/* No iterations were performed. */
/* 1 The equations A*x = b are probably */
/* compatible. Norm(A*x - b) is sufficiently */
/* small, given the values of ATOL and BTOL. */
/* 2 The system A*x = b is probably not */
/* compatible. A least-squares solution has */
/* been obtained that is sufficiently accurate, */
/* given the value of ATOL. */
/* 3 An estimate of cond(Abar) has exceeded */
/* CONLIM. The system A*x = b appears to be */
/* ill-conditioned. Otherwise, there could be an */
/* error in subroutine APROD. */
/* 4 The equations A*x = b are probably */
/* compatible. Norm(A*x - b) is as small as */
/* seems reasonable on this machine. */
/* 5 The system A*x = b is probably not */
/* compatible. A least-squares solution has */
/* been obtained that is as accurate as seems */
/* reasonable on this machine. */
/* 6 Cond(Abar) seems to be so large that there is */
/* no point in doing further iterations, */
/* given the precision of this machine. */
/* There could be an error in subroutine APROD. */
/* 7 The iteration limit ITNLIM was reached. */
/* ITN output The number of iterations performed. */
/* ANORM output An estimate of the Frobenius norm of Abar. */
/* This is the square-root of the sum of squares */
/* of the elements of Abar. */
/* If DAMP is small and if the columns of A */
/* have all been scaled to have length 1.0, */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -