OPERATS(3LAS)

Next: PRECOND(3LAS) Up: Manual Pages Previous: MLSOLV(3LAS)

Asgn_VV, AddAsgn_VV, SubAsgn_VV, MulAsgn_VS, Add_VV, Add_MM, Sub_VV, Sub_MM, Mul_SV, Mul_SM, Mul_SO, Mul_VV, Mul_MV, Mul_OV, MulInv_MV, Transp_M, Transp_O, Diag_M, Upper_M, Lower_M, l1Norm_V, l2Norm_V, MaxNorm_V, OrthoRightKer_VQ, OrthoLeftKer_VQ -- basic operations of linear algebra

SYNOPSIS

#include <laspack/operats.h>
 
Vector *Asgn_VV(Vector *V1, Vector *V2); 
Vector *AddAsgn_VV(Vector *V1, Vector *V2); 
Vector *SubAsgn_VV(Vector *V1, Vector *V2); 
Vector *MulAsgn_VS(Vector *V, double S); 
Vector *Add_VV(Vector *V1, Vector *V2); 
QMatrix *Add_QQ(QMatrix *Q1, QMatrix *Q2); 
Vector *Sub_VV(Vector *V1, Vector *V2); 
QMatrix *Sub_QQ(QMatrix *Q1, QMatrix *Q2); 
Vector *Mul_SV(double S, Vector *V); 
Matrix *Mul_SM(double S, Matrix *M); 
QMatrix *Mul_SQ(double S, QMatrix *Q); 
double Mul_VV(Vector *V1, Vector *V2); 
Vector *Mul_MV(Matrix *M, Vector *V); 
Vector *Mul_QV(QMatrix *Q, Vector *V); 
Vector *MulInv_QV(QMatrix *Q, Vector *V); 
Matrix *Transp_M(Matrix *M); 
QMatrix *Transp_Q(QMatrix *Q); 
QMatrix *Diag_Q(QMatrix *Q); 
QMatrix *Upper_Q(QMatrix *Q); 
QMatrix *Lower_Q(QMatrix *Q); 
double l1Norm_V(Vector *V); 
double l2Norm_V(Vector *V); 
double MaxNorm_V(Vector *V); 
Vector *OrthoRightKer_VQ(Vector *V, QMatrix *Q);
Vector *OrthoLeftKer_VQ(Vector *V, QMatrix *Q);

DESCRIPTION

Functions defined in this module perform the following basic operations of linear algebra:

    Asgn_VV(Vector a, Vector b)     --> Vector   :  a = b
    AddAsgn_VV(Vector a, Vector b)  --> Vector   :  a = a + b
    SubAsgn_VV(Vector a, Vector b)  --> Vector   :  a = a - b
    MulAsgn_VS(Vector a, double s)  --> Vector   :  a = s a
    Add_VV(Vector a, Vector b)      --> Vector   :  a + b
    Add_QQ(QMatrix A, QMatrix B)    --> QMatrix  :  A + B
    Sub_VV(Vector a, Vector b)      --> Vector   :  a - b
    Sub_QQ(QMatrix A, QMatrix B)    --> QMatrix  :  A - B
    Mul_SV(double s, Vector a)      --> Vector   :  s a
    Mul_SM(double s, Matrix P)      --> Matrix   :  s P
    Mul_SQ(double s, QMatrix A)     --> QMatrix  :  s A
    Mul_VV(Vector a, Vector b)      --> double   :  a . b
    Mul_MV(Matrix P, Vector a)      --> Vector   :  P a
    Mul_QV(QMatrix A, Vector a)     --> Vector   :  A a
    MulInv_QV(QMatrix A, Vector a)  --> Vector   :  A^{-1} a
    Transp_M(Matrix P)              --> Matrix   :  P^T
    Transp_Q(QMatrix A)             --> QMatrix  :  A^T
    l1Norm_V(Vector a)              --> double   :  || a ||_1
    l2Norm_V(Vector a)              --> double   :  || a ||_2
    MaxNorm_V(Vector a)             --> double   :  || a ||_max

Here a, b, c are vectors, s a scalar, A, B quadratic matrices, and P a rectangular matrix.

The procedures Diag_Q, Upper_Q and Lower_Q return the diagonal, the strictly upper and the strictly lower triangular part of the matrix Q, respectively.

As usual in C, all above operations produce temporary results which can be processed further. For example by following statements

  Asgn_VV(&a, Add_VV(&b, Mul_SV(s,&c)));
  Asgn_VV(&a, Mul_QV(Transp(&B), &c));

the compound operation a = b + s c and the matrix-vector product a = B c could be performed, respectively. Some sophisticated techniques have been developed in order that this flexible approach do not have essential influence on memory requirement as well as CPU time. Operations with only one parameter of type Vector, Matrix, or QMatrix return objects which share memory extensive arrays of vector components and matrix non-zero elements with the original ones. The modification has an effect only on intern auxiliary variables such as scaling factors for products by scalars, combinations or restrictions of matrices as well as the element ordering for transposition of matrices. This saves memory, avoids unnecessary data transfer between processor and memory, and allows to better exploit some special processors capabilities such as e.g. a multiply-add pipe for compound operations.

The procedures OrthoRightKer_VQ and OrthoLeftKer_VQ carry out the orthogonalization of the vector V on the ``right'' and ``left'' null space of the matrix Q, respectively. The vector V will be as distinguished from other routines modified.

FILES

operats.h ... header file
operats.c ... source file

EXAMPLES

The following example demonstrates how the CG method is by means of above operations implemented in LASPack (simplified for preconditioned systems with non-singular matrices).

Vector *CGIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType PrecondProc, double OmegaPrecond)
{
    int Iter;
    double Alpha, Beta, Rho, RhoOld = 0.0;
    double bNorm;
    size_t Dim;
    Vector r, p, q, z;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", Dim, Normal, True);
    V_Constr(&p, "p", Dim, Normal, True);
    V_Constr(&q, "q", Dim, Normal, True);
    if (PrecondProc != NULL)
        V_Constr(&z, "z", Dim, Normal, True);

    if (LASResult() == LASOK) {
        bNorm = l2Norm_V(b);
        
        Iter = 0;
        /* r = b - A * x(i) */
        if (!IsZero(l1Norm_V(x) / Dim))
            Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        else
            Asgn_VV(&r, b);
        if (PrecondProc != NULL) {
            /* preconditioned CG */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGIterId)
                && Iter < MaxIter) {
                Iter++;
                (*PrecondProc)(A, &z, &r, OmegaPrecond);
                Rho = Mul_VV(&r, &z);
                if (Iter == 1) {
                    Asgn_VV(&p, &z);
           } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
           }
                Asgn_VV(&q, Mul_QV(A, &p));
                Alpha = Rho / Mul_VV(&p, &q);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                RhoOld = Rho;
       }
   } else {
            /* plain CG (z = r) */

            ...

   }
    }
    
    V_Destr(&r);
    V_Destr(&p);
    V_Destr(&q);
    if (PrecondProc != NULL)
        V_Destr(&z);

    Q_Unlock(A);
    V_Unlock(x);
    V_Unlock(b);

    return(x);
}

BUGS

Addition and subtraction of matrices are allowed only for matrices which are derived from the same one. This disables indeed general matrix manipulations but operations of such kind are as a rule not applied in solvers of systems of linear equations.

For example, triangular matrices with weighted diagonal used in the SOR method could be build as follows:

  Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A))