📄 symmlq.cpp
字号:
/* symmlq.f -- translated by f2c (version of 16 May 1991 13:06:06).
You must link the resulting object file with the libraries:
-link <S|C|M|L>f2c.lib (in that order)
*/
#include <math.h>
#include "utils_rqi.h"
#include "aprod.h"
#include "msolve.h"
#include "symmlqblas.h"
#include "f2c.h"
double pow_dd(doublereal *ap, doublereal *bp)
{
return (pow(*ap, *bp));
}
/* Table of constant values */
static integer c__1 = 1;
static doublereal c_b4 = 1.;
static doublereal c_b18 = .33333;
/* Subroutine */ int symmlq_(
integer *n,
doublereal *b, doublereal *r1, doublereal *r2, doublereal *v,
doublereal *w, doublereal *x, doublereal *y, doublereal *work,
logical *checka, logical *goodb, logical *precon,
doublereal *shift,
integer *nout, integer *itnlim,
doublereal *rtol,
integer *istop, integer *itn,
doublereal *anorm, doublereal *acond, doublereal *rnorm, doublereal *ynorm,
doublereal *a, doublereal *vwsqrt, doublereal *orthlist, doublereal *macheps,
doublereal *normxlim,
integer *itnmin)
{
/* System generated locals */
integer i__1;
doublereal d__1, d__2;
/* Builtin functions */
/* Local variables */
static doublereal alfa, diag, dbar, beta, gbar, oldb, epsa;
//extern doublereal ddot_();
static doublereal gmin, gmax, zbar, epsr, epsx, beta1;
//extern doublereal dnrm2_();
static integer i;
static doublereal gamma, s, t, delta, z, denom;
//extern /* Subroutine */ int aprod_();
static doublereal bstep;
//extern /* Subroutine */ int dcopy_();
static doublereal epsln;
//extern /* Subroutine */ int daxpy_();
static doublereal tnorm, cs, ynorm2, sn, cgnorm;
//extern /* Subroutine */ int msolve_();
static doublereal snprod, lqnorm, qrnorm, eps, rhs1, rhs2;
/* ------------------------------------------------------------------
*/
/* SYMMLQ is designed to solve the system of linear equations */
/* Ax = b */
/* where A is an n by n symmetric matrix and b is a given vector. */
/* The matrix A is not required to be positive definite. */
/* (If A is known to be definite, the method of conjugate gradients */
/* might be preferred, since it will require about the same number of
*/
/* iterations as SYMMLQ but slightly less work per iteration.) */
/* The matrix A is intended to be large and sparse. It is accessed */
/* by means of a subroutine call of the form */
/* old call aprod ( n, x, y ) */
/* new: call aprod ( n, x, y, A, vwsqrt, work, orthlist ) -rwl
*/
/* which must return the product y = Ax for any given vector x. */
/* More generally, SYMMLQ is designed to solve the system */
/* (A - shift*I) x = b */
/* where shift is a specified scalar value. If shift and b */
/* are suitably chosen, the computed vector x may approximate an */
/* (unnormalized) eigenvector of A, as in the methods of */
/* inverse iteration and/or Rayleigh-quotient iteration. */
/* Again, the matrix (A - shift*I) need not be positive definite. */
/* The work per iteration is very slightly less if shift = 0. */
/* A further option is that of preconditioning, which may reduce */
/* the number of iterations required. If M = C C' is a positive */
/* definite matrix that is known to approximate (A - shift*I) */
/* in some sense, and if systems of the form My = x can be */
/* solved efficiently, the parameters precon and msolve may be */
/* used (see below). When precon = .true., SYMMLQ will */
/* implicitly solve the system of equations */
/* P (A - shift*I) P' xbar = P b, */
/* i.e. Abar xbar = bbar */
/* where P = C**(-1), */
/* Abar = P (A - shift*I) P', */
/* bbar = P b, */
/* and return the solution x = P' xbar. */
/* The associated residual is rbar = bbar - Abar xbar */
/* = P (b - (A - shift*I)x) */
/* = P r. */
/* In the discussion below, eps refers to the machine precision. */
/* eps is computed by SYMMLQ. A typical value is eps = 2.22e-16 */
/* for IBM mainframes and IEEE double-precision arithmetic. */
/* Parameters */
/* ---------- */
/* n input The dimension of the matrix A. */
/* b(n) input The rhs vector b. */
/* r1(n) workspace */
/* r2(n) workspace */
/* v(n) workspace */
/* w(n) workspace */
/* x(n) output Returns the computed solution x. */
/* y(n) workspace */
/* aprod external A subroutine defining the matrix A. */
/* For a given vector x, the statement */
/* old: call aprod ( n, x, y, ) */
/* new: call aprod ( n, x, y, A, vwsqrt, work, orthl
ist ) -rwl*/
/* must return the product y = Ax */
/* without altering the vector x. */
/* msolve external An optional subroutine defining a */
/* preconditioning matrix M, which should */
/* approximate (A - shift*I) in some sense. */
/* M must be positive definite. */
/* For a given vector x, the statement */
/* old: call msolve( n, x, y ) */
/* new: call msolve( n, x, y ) -rwl */
/* must solve the linear system My = x */
/* without altering the vector x. */
/* In general, M should be chosen so that Abar has
*/
/* clustered eigenvalues. For example, */
/* if A is positive definite, Abar would ideally */
/* be close to a multiple of I. */
/* If A or A - shift*I is indefinite, Abar might */
/* be close to a multiple of diag( I -I ). */
/* NOTE. The program calling SYMMLQ must declare
*/
/* aprod and msolve to be external. */
/* checka input If checka = .true., an extra call of aprod will
*/
/* be used to check if A is symmetric. Also, */
/* if precon = .true., an extra call of msolve */
/* will be used to check if M is symmetric. */
/* goodb input Usually, goodb should be .false. */
/* If x is expected to contain a large multiple of
*/
/* b (as in Rayleigh-quotient iteration), */
/* better precision may result if goodb = .true. */
/* See Lewis (1977) below. */
/* When goodb = .true., an extra call to msolve */
/* is required. */
/* precon input If precon = .true., preconditioning will */
/* be invoked. Otherwise, subroutine msolve */
/* will not be referenced; in this case the */
/* actual parameter corresponding to msolve may */
/* be the same as that corresponding to aprod. */
/* shift input Should be zero if the system Ax = b is to be */
/* solved. Otherwise, it could be an */
/* approximation to an eigenvalue of A, such as */
/* the Rayleigh quotient b'Ab / (b'b) */
/* corresponding to the vector b. */
/* If b is sufficiently like an eigenvector */
/* corresponding to an eigenvalue near shift, */
/* then the computed x may have very large */
/* components. When normalized, x may be */
/* closer to an eigenvector than b. */
/* nout input A file number. */
/* If nout .gt. 0, a summary of the iterations */
/* will be printed on unit nout. */
/* itnlim input An upper limit on the number of iterations. */
/* rtol input A user-specified tolerance. SYMMLQ terminates
*/
/* if it appears that norm(rbar) is smaller than */
/* rtol * norm(Abar) * norm(xbar), */
/* where rbar is the transformed residual vector,
*/
/* rbar = bbar - Abar xbar. */
/* If shift = 0 and precon = .false., SYMMLQ */
/* terminates if norm(b - A*x) is smaller than */
/* rtol * norm(A) * norm(x). */
/* istop output An integer giving the reason for termination...
*/
/* -1 beta2 = 0 in the Lanczos iteration; i.e. the */
/* second Lanczos vector is zero. This means the
*/
/* rhs is very special. */
/* If there is no preconditioner, b is an */
/* eigenvector of A. */
/* Otherwise (if precon is true), let My = b. */
/* If shift is zero, y is a solution of the */
/* generalized eigenvalue problem Ay = lambda My,
*/
/* with lambda = alpha1 from the Lanczos vectors.
*/
/* In general, (A - shift*I)x = b */
/* has the solution x = (1/alpha1) y */
/* where My = b. */
/* 0 b = 0, so the exact solution is x = 0. */
/* No iterations were performed. */
/* 1 Norm(rbar) appears to be less than */
/* the value rtol * norm(Abar) * norm(xbar). */
/* The solution in x should be acceptable. */
/* 2 Norm(rbar) appears to be less than */
/* the value eps * norm(Abar) * norm(xbar). */
/* This means that the residual is as small as */
/* seems reasonable on this machine. */
/* 3 Norm(Abar) * norm(xbar) exceeds norm(b)/eps, */
/* which should indicate that x has essentially */
/* converged to an eigenvector of A */
/* corresponding to the eigenvalue shift. */
/* 4 acond (see below) has exceeded 0.1/eps, so */
/* the matrix Abar must be very ill-conditioned. */
/* x may not contain an acceptable solution. */
/* 5 The iteration limit was reached before any of */
/* the previous criteria were satisfied. */
/* 6 The matrix defined by aprod does not appear */
/* to be symmetric. */
/* For certain vectors y = Av and r = Ay, the */
/* products y'y and r'v differ significantly. */
/* 7 The matrix defined by msolve does not appear */
/* to be symmetric. */
/* For vectors satisfying My = v and Mr = y, the */
/* products y'y and r'v differ significantly. */
/* 8 An inner product of the form x' M**(-1) x */
/* was not positive, so the preconditioning matrix
*/
/* M does not appear to be positive definite. */
/* If istop .ge. 5, the final x may not be an */
/* acceptable solution. */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -