Digital

Name

ditsol_plscg − Preconditioned least square conjugate gradient method (Serial and Parallel Versions)

FORMAT

DITSOL_PLSCG (matvec, pcondl, pcondr, mstop, a, ia, x, b, n,
              ql, iql , qr, iqr, iparam, rparam, iwork, rwork, ierror)

Arguments

DITSOL_PLSCG has the standard parameter list for an iterative solver. 

Description

The least squares conjugate gradient is a robust method for the  solution of general linear systems. It is equivalent to applying the conjugate gradient method to the normal equations:

     transp(A)∗ A ∗ x = transp(A)∗ b

This  method requires the evaluation of two matrix products, involving matrix A and transp(A). It suffers from the drawback that the condition number of transp(A) ∗ A is the square of the condition number of A, and therefore the convergence of the method is slow. To alleviate the numerical instability  resulting from a straightforward application of the conjugate  gradient method to the normal equations, DXML adopts the  implementation proposed in [Bjorck and Elfving 1979]. 

The implementation of the least squares conjugate gradient  method requires the routine MATVEC to provide operations for  both job= 0 and job= 1. The routines MATVEC, PCONDL (if used), PCONDR (if used) and MSTOP (if used) should be declared external in your calling (sub)program. 

DXML provides the following four forms of the method:

•Unpreconditioned least squares conjugate gradient method:

This is the conjugate gradient method applied to

    transp(A) ∗ A ∗ x = transp(A) ∗ b

where A is a general matrix. As no preconditioning is used,  both PCONDL and PCONDR are dummy input parameters. 

For the unpreconditioned least squares conjugate gradient  method, the length of the real work space array, defined by the variable nrwk (IPARAM(4)), should be at least 4∗n, where  n is the order of the matrix A. 

The vector z, passed as an input argument to the routine  MSTOP, is not defined. 

•Least squares conjugate gradient method with left preconditioning:

This is the conjugate gradient method applied to

   ( transp(A) ∗ inverse(transp(QL)) ∗  inverse(QL)  ∗ A) ∗ x  =
       (transp(A) ∗ inverse(transp(QL)) ∗ inverse(QL) ∗ b )

The routine PCONDL, with job= 0 should evaluate

   v = inverse(QL) ∗ u

and with job= 1 should evaluate

   v = inv_transp(QL) ∗ u

The routine PCONDR is not used and is therefore a dummy input parameter. 

For the least squares conjugate gradient method, with left preconditioning, the length of the real work space array,  defined by the variable nrwk (IPARAM(4)), should be at  least 5∗n, where n is the order of the matrix A. This does  not include the memory requirements of the preconditioner. 

The vector z, passed as an input argument to the routine  MSTOP, is defined as

    z = inverse(QL) ∗ r

where r is the residual at the i-th iteration. 

•Least squares conjugate gradient method with right preconditioning:

This is the conjugate gradient method applied to

   ( inverse(transp(QR)) ∗ transp(A) ∗ A ∗ inverse(QR) ) ∗ y  =
         ( inverse(transp(QR)) ∗ transp(A) ∗ b )

where

   y = QR ∗ x

The routine PCONDR, with job= 0 should evaluate

   v = inverse(QR) ∗ u

and with job= 1 should evaluate

   v = inv_transp(QR) ∗ u

The routine PCONDL is not used and is therefore a dummy input parameter. 

For the least squares conjugate gradient method, with right  preconditioning, the length of the real work space array,  defined by the variable nrwk (IPARAM(4)), should be at least  5∗n, where n is the order of the matrix A. This does not  include the memory requirements of the preconditioner. 

The vector z, passed as an input argument to the routine  MSTOP, is not defined. 

•Least squares conjugate gradient method with split preconditioning:

This is the conjugate gradient method applied to

   ( inverse(transp(QR)) ∗ transp(A) ∗ inverse(transp(QL) ∗
       inverse(QL)  ∗ A ∗ inverse(QR)) ∗ y  =
       (inverse(transp(QR)) ∗ transp(A) ∗ inverse(transp(QL)) ∗
       inverse(QL) ∗ b)

where

   y = QR ∗ x

The routine PCONDL, with job= 0 should evaluate

   v = inverse(QL) ∗ u

and with job= 1 should evaluate

   v = inv_transp(QL)∗ u

The routine PCONDR, with job= 0 should evaluate

   v = inverse(QR) ∗ u

and with job= 1 should evaluate

   v = inv_transp(QR) ∗ u

For the least squares conjugate gradient method, with split  preconditioning, the length of the real work space array, defined by the variable nrwk (IPARAM(4)), should be at least  6∗n, where n is the order of the matrix A. This does not  include the memory requirements of the preconditioner. 

The vector z, passed as an input argument to the routine MSTOP, is defined as

   z = inverse(QL) ∗ r

where r is the residual at the i-th iteration. 

This routine is available in both serial and parallel versions. The routine names and parameter list are identical for both versions. For information about linking to the serial or to the parallel library, refer to the DXML Reference Manual.