guslv.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 guslv                                                          */
/*****************************************************************************/
/*****************************************************************************/
/*******     guslv ..... GENERAL UPPER 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 guslv(int *neqns, int *xlnz, double *lnz, int *xnzsub, int *nzsub, 
	   double *diag, double *rhs)
{
   int    k,*pn;
   double *p,*pl,*pr;

   /* ------------------------- */
   /* BACKWARD SUBSTITUTION ... */
   /* ------------------------- */
   /* use STACK VARIABLES xlnz,xnzsub,diag as pointers. The real variables   */
   /* are NOT AFFECTED by these changes!                                     */
   xlnz  +=*neqns;
   xnzsub+=*neqns;
   diag  +=*neqns;
   /* adjust rhs,lnz, nzsub */
   rhs--; lnz--; nzsub--;
   pr=rhs+*neqns;
   while (pr>rhs)
   {
         xlnz--; xnzsub--; diag--;
	 pn=nzsub+*xnzsub;
	 pl=lnz+xlnz[0];
	 k=(xlnz[1]-xlnz[0])%4;
	 p =lnz+xlnz[1]-k;
	 for (; pl<p; pl+=4,pn+=4)
	 {
	     *pr-=pl[0]*rhs[pn[0]]+pl[1]*rhs[pn[1]]
	         +pl[2]*rhs[pn[2]]+pl[3]*rhs[pn[3]];
	 } /* end for */
	 p+=k;
	 while (pl<p)
	       *pr-=*pl++*rhs[*pn++];
	 *pr--/=*diag;
   } /* end for */
   /* not necessary: re-adjust rhs,lnz, nzsub */
} /* end guslv */