glslv.c

http://www.math.tu-berlin.de/~bolle/SPARSPAK/index-e.html

#include <sparspak.h>
/* #include <stdlib.h> */
/* #include <string.h> */
/* #define PRINT_INFO */

/* #define MC_EXTERNAL       0 */
/* #define Malloc(A,B)       malloc(A) */
/* #define Realloc(A,B,C)    realloc(A,B) */
/* #define Free(A)           free(A) */


/* ----- void glslv                                                          */
/*****************************************************************************/
/*****************************************************************************/
/*******     glslv ..... GENERAL LOWER TRIANGULAR SPARSE SOLVE    ************/
/*****************************************************************************/
/*****************************************************************************/
/*    PURPOSE - TO PERFORM SOLUTION OF A FACTORED SYSTEM, WHERE              */
/*       THE MATRIX IS STORED IN THE COMPRESSED SUBSCRIPT                    */
/*       SPARSE FORMAT.                                                      */
/*                                                                           */
/*    INPUT PARAMETERS -                                                     */
/*       neqns           - NUMBER OF EQUATIONS.                              */
/*       (xlnz, lnz)     - STRUCTURE OF NONZEROS IN L.                       */
/*       (xnzsub, nzsub) - COMPRESSED SUBSCRIPT STRUCTURE.                   */
/*       diag            - DIAGONAL COMPONENTS OF L                          */
/*                                                                           */
/*    UPDATED PARAMETER -                                                    */
/*       rhs - ON INPUT, IT CONTAINS THE RHS VECTOR, AND ON                  */
/*              OUTPUT, THE SOLUTION VECTOR.                                 */
/*                                                                           */
/****************************************************************            */
/*                                                                           */
void glslv(int *neqns, int *xlnz, double *lnz, int *xnzsub, int *nzsub, 
	    double *diag, double *rhs)
{
   int     k, *pn;
   double  *pr,*p;

   /* ------------------------ */
   /* FORWARD SUBSTITUTION ... */
   /* ------------------------ */
   /* use STACK VARIABLES 'lnz,xlnz,xnzsub,diag' as pointers */
   lnz+=xlnz[0]-1;
   pr=rhs;
   rhs--; nzsub--;
   for (; pr<=rhs+*neqns; xlnz++,xnzsub++,diag++,pr++)
   {
       *pr/=*diag;
       pn=nzsub+*xnzsub;
       k=(xlnz[1]-xlnz[0])%4;
       p=lnz+(xlnz[1]-xlnz[0])-k;
       for (; lnz<p; lnz+=4,pn+=4)
       {
	   rhs[pn[0]]-=lnz[0]**pr;   rhs[pn[1]]-=lnz[1]**pr;
	   rhs[pn[2]]-=lnz[2]**pr;   rhs[pn[3]]-=lnz[3]**pr;
       } /* end for */
       p+=k;
       while (lnz<p)
	     rhs[*pn++]-=*lnz++**pr;
   } /* end for */
} /* end glslv */