next up previous contents index
Next: PRECOND(3LAS) Up: Manual Pages Previous: MLSOLV(3LAS)

  OPERATS(3LAS)

NAME

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 is a scalar, A, B are quadratic matrices, and P is 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 the above operations produce temporary results which can be processed further. For example by the 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 , respectively, could be performed. Some sophisticated techniques have been developed in order that this flexible approach may not have an essential influence on the memory requirement as well as the CPU time. Operations with only one parameter of the 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 internal 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 us to better exploit some special processor 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 modified as distinguished from other routines.

FILES

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

EXAMPLES

The following example demonstrates how the CG method is by means of the 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);
}

SEE ALSO

vector(3LAS), matrix(3LAS), qmatrix(3LAS), errhandl(3LAS)

BUGS

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

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

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



next up previous contents index
Next: PRECOND(3LAS) Up: Manual Pages Previous: MLSOLV(3LAS)



Tomas Skalicky (skalicky@msmfs1.mw.tu-dresden.de)