chlsetup.c

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

#include <stdio.h>

#include <stdlib.h>
#include <string.h>
#include <hbo.h>
#include <sparspak.h>
#define PRINT_INFO

/* void chlsetup                                                             */
/* PURPOSE                                                                   */
/*        SETUP DATA FOR THE CHOLESKY FACTOR. THE CHOLESKY FACTOR IS         */
/*        INITIALIZED WITH THE DATA OF A. ZEROS FOR LATER FILL-IN ARE        */
/*        ALREADY INCLUDED.                                                  */
/*        FROM 'smbfct' WE KNOW, WHERE WE HAVE TO INSERT ZEROS FOR LATER     */
/*        FILL-IN. SO IT IS NECESSARY TO RUN 'smbfct' FIRST.                 */
/*                                                                           */
/*     INPUT PARAMETERS -                                                    */
/*        neqns           - NUMBER OF EQUATIONS.                             */
/*        xlnz            - INDEX INTO THE NONZERO STORAGE VECTOR lnz.       */
/*        (xnzsub, nzsub) - THE COMPRESSED SUBSCRIPT VECTORS.                */
/*        (matptr, matx)  - SUBSCRIPT VECTORS OF A                           */
/*        matval          - DATA OF A                                        */
/*        (xadj, adjncy)  - THE ADJACENCY STRUCTURE.                         */
/*        (perm, invp)    - THE PERMUTATION VECTOR AND ITS INVERSE.          */
/*                                                                           */
/*     OUTPUT PARAMETERS -                                                   */
/*        lnz             - STORAGE VECTOR                                   */
/*        diag            - THE DIAGONAL OF A.                               */
void chlsetup(int *neqns, int *xlnz, double *lnz, int *xnzsub, int *nzsub,
	      double *diag, int *matptr, int *matx, double *matval, 
	      int *xadj, int *adjncy, int *perm, int *invp)
{
   int    flag,
          i, j, j2, j3, k;

   /* extract diagonal entry */
   for (i=0; i<*neqns; i++)
       diag[i]=matval[matptr[perm[i]-1]-1]; 
   
   k=0; 
   for (i=0; i<*neqns-1; i++) 
   { 
       j2=xnzsub[i]-1;
       while (j2<xnzsub[i]-1+xlnz[i+1]-xlnz[i]) 
       {
	     /* compare new entry with adjacency list of node perm[i] */
	     flag=0;
	     j=xadj[perm[i]-1]-1;
	     while (j<xadj[perm[i]]-1 && flag==0)
	     { 
		   if (invp[adjncy[j]-1]==nzsub[j2])
		      flag=-1;
		   else
		      j++; 		
	     } /* end while */ 

	     /* copy old entry */
	     if (flag)
	     {
	        /* search corresponding value in 'matx' */
	        /* search in column perm[i] */
	        if (adjncy[j]>perm[i])
		{
		   /* start with first off-diagonal entry */
		   j3=matptr[perm[i]-1];
		   while (matx[j3]<adjncy[j])
		         j3++;
		} /* end if */
		/* search in column adjncy[j] */
		else
		{
		   /* start with first off-diagonal entry */
		   j3=matptr[adjncy[j]-1];
		   while (matx[j3]<perm[i])
		         j3++; 
		} /* end else */
		lnz[k++]=matval[j3];
	     } /* end if */
	     /* insert zero for later fill-in */
	     else
	        lnz[k++]=0;
	     j2++;
       } /* end while */ 
   } /* end for */ 
}