/* * TEST MODULE for the sparse matrix routines * * Author: Advisor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This file contains the test routine for the sparse matrix routines. * They are able to read matrices from files and solve them. * * >>> Functions contained in this file: * main * ReadMatrixFromFile * PrintMatrixErrorMessage * InterpretCommandLine * GetProgramName * Usage */ /* * Revision and copyright information. * * Copyright (c) 1985,86,87,88 * by Kenneth S. Kundert and the University of California. * * Permission to use, copy, modify, and distribute this software and * its documentation for any purpose and without fee is hereby granted, * provided that the copyright notices appear in all copies and * supporting documentation and that the authors and the University of * California are properly credited. The authors and the University of * California make no representations as to the suitability of this * software for any purpose. It is provided `as is', without express * or implied warranty. */ #ifndef lint static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; static char RCSid[] = "@(#)$Header: spTest.c,v 1.2 88/06/18 11:16:24 kundert Exp $"; #endif /* * IMPORTS * * >>> Import descriptions: * stdio.h math.h ctype.h * Standard C libraries. * spConfig.h * Macros that customize the sparse matrix package. It is not normally * necessary, nor is normally particularly desirable to include this * file into the calling routines. Nor should spINSIDE_SPARSE be defined. * It is done in this test file so that the complex test routines may be * removed when they don't exist in Sparse. * spmatrix.h * Macros and declarations to be imported by the user. */ #define YES 1 #define NO 0 #include #include #include #define spINSIDE_SPARSE #include "spConfig.h" #undef spINSIDE_SPARSE #include "spmatrix.h" /* Declarations that should be in stdio.h. */ extern char *malloc(); #ifdef ultrix extern void free(); extern void exit(); extern void perror(); #else extern free(); extern exit(); extern perror(); #endif /* Declarations that should be in strings.h. */ extern int strcmp(), strncmp(), strlen(); #ifdef lint #undef MULTIPLICATION #undef DETERMINANT #define MULTIPLICATION YES #define DETERMINANT YES #define LINT YES #else /* not lint */ #define LINT NO #endif /* not lint */ /* * TIMING */ #ifndef vms # include #endif #ifdef notdef /* __STDC__ */ /* Deleted because some ANSI C compilers to not yet have complete h files. */ # include # define HZ CLK_TCK #endif #ifndef HZ # ifdef vax # ifdef unix # define HZ 60 # endif # ifdef vms # define HZ 100 # endif # endif # ifdef hpux # ifdef hp300 # define HZ 50 # endif # ifdef hp500 # define HZ 60 # endif # endif #endif /* Routine that queries the system to find the process time. */ double Time() { struct time {long user, system, childuser, childsystem;} time; (void)times(&time); return (double)time.user / (double)HZ; } /* * MACROS */ #define BOOLEAN int #define NOT ! #define AND && #define OR || #define ALLOC(type,number) ((type *)malloc((unsigned)(sizeof(type)*(number)))) #define ABS(a) ((a) < 0.0 ? -(a) : (a)) #ifdef MAX #undef MAX #endif #ifdef MIN #undef MIN #endif #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MIN(a,b) ((a) < (b) ? (a) : (b)) /* * IMAGINARY VECTORS * * The imaginary vectors iRHS and iSolution are only needed when the * options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set. The following * macro makes it easy to include or exclude these vectors as needed. */ #if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS #define IMAG_VECTORS , iRHS, iSolution #else #define IMAG_VECTORS #endif /* * GLOBAL DECLARATIONS * * These variables, types, and macros are used throughout this file. * * >>> Macros * PRINT_LIMIT * The maximum number of terms to be printed form the solution vector. * May be overridden by the -n option. * * >>> Variable types: * ElementRecord * A structure for holding data for each matrix element. * * >>> Global variables: * ProgramName (char *) * The name of the test program with any path stripped off. * FileName (char *) * The name of the current file name. * Complex (BOOLEAN) * The flag that indicates whether the matrix is complex or real. * Element (ElementRecord[]) * Array used to hold information about all the elements. * Matrix (char *) * The pointer to the matrix. * Size (int) * The size of the matrix. * RHS (RealVector) * The right-hand side vector (b in Ax = b). * Solution (RealVector) * The solution vector (x in Ax = b). * RHS_Verif (RealVector) * The calculated RHS vector. */ /* Begin global declarations. */ #define PRINT_LIMIT 9 typedef spREAL RealNumber, *RealVector; typedef struct { RealNumber Real; RealNumber Imag; } ComplexNumber, *ComplexVector; static char *ProgramName; static char *FileName; static int Size, MaxSize = 0; static int PrintLimit = PRINT_LIMIT, Iterations = 1, ColumnAsRHS = 1; static BOOLEAN Complex, SolutionOnly = NO, RealAsComplex = NO, Transposed = NO; static BOOLEAN CreatePlotFiles = NO, UseColumnAsRHS = NO, PrintLimitSet = NO; static BOOLEAN ExpansionWarningGiven, DiagPivoting = DIAGONAL_PIVOTING; static char Message[BUFSIZ], *Matrix = NULL; static RealNumber AbsThreshold = (RealNumber)0, RelThreshold = (RealNumber)0; static RealVector RHS, Solution, RHS_Verif; static RealVector iRHS, iSolution, iRHS_Verif; static ComplexVector cRHS, cSolution, cRHS_Verif; int Init(); /* * MAIN TEST ROUTINE for sparse matrix routines * * This routine reads a matrix from a file and solves it several * times. The solution is printed along with some statistics. * The format expected for the input matrix is the same as what is output by * spFileMatrix and spFileVector. * * >>> Local variables: * Determinant (RealNumber) * Real portion of the determinant of the matrix. * iDeterminant (RealNumber) * Imaginary portion of the determinant of the matrix. * Error (int) * Holds error status. * Last (int) * The index of the last term to be printed of the solution. * Iterations (int) * The number of times that the matrix will be factored and solved. * MaxRHS (RealNumber) * The largest term in the given RHS vector. * Residual (RealNumber) * The sum of the magnitude of the differences in each corresponding * term of the given and calculated RHS vector. * Exponent (int) * Exponent for the determinant. * Growth (RealNumber) * The growth that has occurred in the matrix during the factorization. * * >>> Timing variables: * BuildTime (double) * The time required to build up the matrix including the time to clear * the matrix. * ConditionTime (double) * The time required to compute the condition number. * DeterminantTime (double) * The time required to compute the determinant. * FactorTime (double) * The time required to factor the matrix without ordering. * InitialFactorTime (double) * The time required to factor the matrix with ordering. * SolveTime (double) * The time required to perform the forward and backward elimination. * StartTime (double) * The time that a timing interval was started. * * >>> Global variables: * Complex * Matrix * Size * RHS * iRHS */ main(ArgCount, Args) int ArgCount; char *Args[]; { register long I; int Error, Last, Elements, Fillins; RealNumber MaxRHS, Residual; RealNumber Determinant, iDeterminant; int Exponent; RealNumber ConditionNumber, PseudoCondition; RealNumber LargestBefore, LargestAfter, Roundoff, InfNorm; BOOLEAN StandardInput; char j, PlotFile[BUFSIZ], ErrMsg[BUFSIZ]; double StartTime, BeginTime, BuildTime, FactorTime, SolveTime, PartitionTime; double InitialFactorTime, ConditionTime, DeterminantTime; extern double Time(); extern char *sbrk(); /* Begin `main'. */ BeginTime = Time(); ArgCount = InterpretCommandLine( ArgCount, Args ); /* Assure that the Sparse is compatible with this test program.*/ # if NOT EXPANDABLE OR NOT INITIALIZE OR NOT ARRAY_OFFSET fprintf(stderr, "%s: Sparse is configured inappropriately for test program.\n", ProgramName); fprintf(stderr, " Enable EXPANDABLE, INITIALIZE, and ARRAY_OFFSET in `spConfig.h'.\n"); exit(1); # endif do { /* Initialization. */ BuildTime = FactorTime = SolveTime = 0.0; ExpansionWarningGiven = NO; /* Create matrix. */ Matrix = spCreate( 0, spCOMPLEX, &Error ); if (Matrix == NULL) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } if( Error >= spFATAL ) goto End; /* Read matrix. */ if (ArgCount == 0) FileName = NULL; else FileName = *(++Args); Error = ReadMatrixFromFile(); if( Error) goto End; StandardInput = (FileName == NULL); /* Clear solution vector if row and column numbers are not densely packed. */ if (spGetSize(Matrix, YES) != spGetSize(Matrix, NO)) { if (Complex OR RealAsComplex) { for (I = Size; I > 0; I--) cSolution[I].Real = cSolution[I].Imag = 0.0; } else { for (I = Size; I > 0; I--) Solution[I] = 0.0; } } /* Perform initial build, factor, and solve. */ (void)spInitialize(Matrix, Init); #if MODIFIED_NODAL spMNA_Preorder( Matrix ); #endif #if DOCUMENTATION if (CreatePlotFiles) { if (StandardInput) (void)sprintf(PlotFile, "bef"); else (void)sprintf(PlotFile, "%s.bef", FileName); if (NOT spFileMatrix( Matrix, PlotFile, FileName, NO, NO, NO )) { (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile); perror( ErrMsg ); } } #if NO spPrint( Matrix, NO /*reodered*/, NO /*data*/, NO /*header*/ ); #endif #endif /* DOCUMENTATION */ #if STABILITY if (NOT SolutionOnly) LargestBefore = spLargestElement(Matrix); #endif #if CONDITION if (NOT SolutionOnly) InfNorm = spNorm(Matrix); #endif StartTime = Time(); Error = spOrderAndFactor( Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting ); InitialFactorTime = Time() - StartTime; PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) goto End; #if DOCUMENTATION if (CreatePlotFiles) { if (StandardInput) (void)sprintf(PlotFile, "aft"); else (void)sprintf(PlotFile, "%s.aft", FileName); if (NOT spFileMatrix( Matrix, PlotFile, FileName, YES, NO, NO )) { (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile); perror( ErrMsg ); } } #if NO spFileStats( Matrix, FileName, "stats" ); #endif #endif /* DOCUMENTATION */ /* * IMAG_VECTORS is a macro that replaces itself with `, iRHS, iSolution' * if the options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set, * otherwise it disappears without a trace. */ #if TRANSPOSE if (Transposed) spSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS ); else #endif spSolve( Matrix, RHS, Solution IMAG_VECTORS ); if (SolutionOnly) Iterations = 0; else { #if STABILITY LargestAfter = spLargestElement(Matrix); Roundoff = spRoundoff(Matrix, LargestAfter); #endif #if CONDITION StartTime = Time(); ConditionNumber = spCondition(Matrix, InfNorm, &Error); ConditionTime = Time() - StartTime; PrintMatrixErrorMessage( Error ); #endif #if PSEUDOCONDITION PseudoCondition = spPseudoCondition(Matrix); #endif StartTime = Time(); spPartition( Matrix, spDEFAULT_PARTITION ); PartitionTime = Time() - StartTime; } /* Solve system of equations Iterations times. */ for(I = 1; I <= Iterations; I++) { StartTime = Time(); (void)spInitialize(Matrix, Init); BuildTime += Time() - StartTime; StartTime = Time(); Error = spFactor( Matrix ); FactorTime += Time() - StartTime; if( Error != spOKAY ) PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) goto End; StartTime = Time(); /* * IMAG_VECTORS is a macro that replaces itself with `, iRHS, iSolution' * if the options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set, * otherwise it disappears without a trace. */ #if TRANSPOSE if (Transposed) spSolveTransposed(Matrix, RHS, Solution IMAG_VECTORS ); else #endif spSolve( Matrix, RHS, Solution IMAG_VECTORS ); SolveTime += Time() - StartTime; } /* Print Solution. */ if (SolutionOnly) { if (PrintLimitSet) Last = MIN( PrintLimit, Size ); else Last = Size; j = ' '; } else { Last = (PrintLimit > Size) ? Size : PrintLimit; if (Last > 0) printf("Solution:\n"); j = 'j'; } if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Last; I++) { printf("%-16.9lg %-.9lg%c\n", (double)Solution[I], (double)iSolution[I], j); } #else for (I = 1; I <= Last; I++) { printf("%-16.9lg %-.9lg%c\n", (double)cSolution[I].Real, (double)cSolution[I].Imag, j); } #endif } else { for (I = 1; I <= Last; I++) printf("%-.9lg\n", (double)Solution[I]); } if (Last < Size AND Last != 0) printf("Solution list truncated.\n"); printf("\n"); #if DETERMINANT /* Calculate determinant. */ if (NOT SolutionOnly) { StartTime = Time(); #if spCOMPLEX spDeterminant( Matrix, &Exponent, &Determinant, &iDeterminant ); #else spDeterminant( Matrix, &Exponent, &Determinant ); #endif DeterminantTime = Time() - StartTime; if (Complex OR RealAsComplex) { Determinant = hypot( Determinant, iDeterminant ); while (Determinant >= 10.0) { Determinant *= 0.1; Exponent++; } } } #else Determinant = 0.0; Exponent = 0; #endif #if MULTIPLICATION if (NOT SolutionOnly) { /* Calculate difference between actual RHS vector and RHS vector * calculated from solution. */ /* Find the largest element in the given RHS vector. */ MaxRHS = 0.0; if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) { if (ABS(RHS[I]) > MaxRHS) MaxRHS = ABS(RHS[I]); if (ABS(iRHS[I]) > MaxRHS) MaxRHS = ABS(iRHS[I]); } #else for (I = 1; I <= Size; I++) { if (ABS(cRHS[I].Real) > MaxRHS) MaxRHS = ABS(cRHS[I].Real); if (ABS(cRHS[I].Imag) > MaxRHS) MaxRHS = ABS(cRHS[I].Imag); } #endif } else { for (I = 1; I <= Size; I++) { if (ABS(RHS[I]) > MaxRHS) MaxRHS = ABS(RHS[I]); } } /* Rebuild matrix. */ (void)spInitialize(Matrix, Init); #if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS if (Transposed) { spMultTransposed( Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution); } else spMultiply(Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution); #else if (Transposed) spMultTransposed( Matrix, RHS_Verif, Solution ); else spMultiply( Matrix, RHS_Verif, Solution ); #endif /* Calculate residual. */ Residual = 0.0; if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) { Residual += ABS(RHS[I] - RHS_Verif[I]) + ABS(iRHS[I] - iRHS_Verif[I]); } #else for (I = 1; I <= Size; I++) { Residual += ABS(cRHS[I].Real - cRHS_Verif[I].Real) + ABS(cRHS[I].Imag - cRHS_Verif[I].Imag); } #endif } else { for (I = 1; I <= Size; I++) Residual += ABS(RHS[I] - RHS_Verif[I]); } } #endif /* Print statistics. */ if (NOT SolutionOnly) { Elements = spElementCount(Matrix); Fillins = spFillinCount(Matrix); printf("Initial factor time = %.2lf.\n", InitialFactorTime); printf("Partition time = %.2lf.\n", PartitionTime); if (Iterations > 0) { printf("Build time = %.3lf.\n", (BuildTime/Iterations)); printf("Factor time = %.3lf.\n",(FactorTime/Iterations)); printf("Solve time = %.3lf.\n", (SolveTime/Iterations)); } #if STABILITY printf("Condition time = %.2lf.\n", ConditionTime); #endif #if DETERMINANT printf("Determinant time = %.2lf.\n", DeterminantTime); #endif printf("\nTotal number of elements = %d.\n", Elements); printf("Average number of elements per row initially = %.2lf.\n", (double)(Elements - Fillins) / (double)spGetSize(Matrix, NO)); printf("Total number of fill-ins = %d.\n", Fillins); #if DETERMINANT OR MULTIPLICATION OR PSEUDOCONDITION OR CONDITION OR STABILITY putchar('\n'); #endif #if STABILITY if (LargestBefore != 0.0) printf("Growth = %.2lg.\n", LargestAfter / LargestBefore); printf("Max error in matrix = %.2lg.\n", Roundoff); #endif #if STABILITY if(ABS(ConditionNumber) > 10 * SMALLEST_REAL); printf("Condition number = %.2lg.\n", 1.0 / ConditionNumber); #endif #if CONDITION AND STABILITY printf("Estimated upper bound of error in solution = %.2lg.\n", Roundoff / ConditionNumber); #endif #if PSEUDOCONDITION printf("PseudoCondition = %.2lg.\n", PseudoCondition); #endif #if DETERMINANT printf("Determinant = %.3lg", (double)Determinant ); if (Determinant != 0.0 AND Exponent != 0) printf("e%d", Exponent); putchar('.'); putchar('\n'); #endif #if MULTIPLICATION if (MaxRHS != 0.0) printf("Normalized residual = %.2lg.\n", (Residual / MaxRHS)); #endif } End:; (void)fflush(stdout); CheckInAllComplexArrays(); spDestroy(Matrix); Matrix = NULL; } while( --ArgCount > 0 ); if (NOT SolutionOnly) { printf("\nAggregate resource usage:\n"); printf(" Time required = %.2lf seconds.\n", Time() - BeginTime); printf(" Virtual memory used = %d kBytes.\n\n", ((int)sbrk(0))/1000); } exit (0); } /* * READ MATRIX FROM FILE * * This function reads the input file for the matrix and the RHS vector. * If no RHS vector exists, one is created. If there is an error in the * file, the appropriate error messages are delivered to standard output. * * >>> Returned: * The error status is returned. If no error occurred, a zero is returned. * Otherwise, a one is returned. * * >>> Local variables: * pMatrixFile (FILE *) * The pointer to the file that holds the matrix. * InputString (char []) * String variable for holding input from the matrix file. * Message (char []) * String variable that contains a one line descriptor of the matrix. * Descriptor is taken from matrix file. * * >>> Global variables: * Complex * Size * Element * RHS * iRHS */ static int ReadMatrixFromFile() { long I, Reads; FILE *pMatrixFile; char sInput[BUFSIZ], sType[BUFSIZ], *p; int Error, Row, Col, Count = 0, LineNumber, RHS_Col, IntSize; double Real, Imag = 0.0; ComplexNumber *pValue, *pInitInfo, *CheckOutComplexArray(); RealNumber *pElement; static char *EndOfFile = "%s: unexpected end of file `%s' at line %d.\n"; static char *Syntax = "%s: syntax error in file `%s' at line %d.\n"; /* Begin `ReadMatrixFromFile'. */ /* Open matrix file in read mode. */ if (NOT SolutionOnly) putchar('\n'); if (FileName == NULL) { FileName = "standard input"; pMatrixFile = stdin; } else { pMatrixFile = fopen(FileName, "r"); if (pMatrixFile == NULL) { fprintf(stderr, "%s: file %s was not found.\n", ProgramName, FileName); return 1; } } Complex = NO; LineNumber = 1; /* Read and print label. */ if (NULL == fgets( Message, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } /* For compatibility with the old file syntax. */ if (!strncmp( Message, "Starting", 8 )) { /* Test for complex matrix. */ if (strncmp( Message, "Starting complex", 15 ) == 0) Complex = YES; LineNumber++; if (NULL == fgets( Message, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } } if (NOT SolutionOnly) printf("%-s\n", Message); /* Read size of matrix and determine type of matrix. */ LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } if ((Reads = sscanf( sInput,"%d %s", &Size, sType )) < 1) { fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1; } if (Reads == 2) { for (p = sType; *p != '\0'; p++) if (isupper(*p)) *p += 'a'-'A'; if (strncmp( sType, "complex", 7 ) == 0) Complex = YES; else if (strncmp( sType, "real", 7 ) == 0) Complex = NO; else { fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1; } } EnlargeVectors( Size, YES ); RHS_Col = MIN( Size, ColumnAsRHS ); #if NOT spCOMPLEX if (Complex) { fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n", ProgramName); fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n"); return 1; } #endif #if NOT REAL if (NOT (Complex OR RealAsComplex)) { fprintf(stderr, "%s: Sparse is not configured to solve real matrices.\n", ProgramName); fprintf(stderr," Enable REAL in `spConfig.h'.\n"); return 1; } #endif /* Read matrix elements. */ do { if (Count == 0) pValue = CheckOutComplexArray( Count = 1000 ); LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } if (Complex) { Reads = sscanf( sInput,"%d%d%lf%lf", &Row, &Col, &Real, &Imag ); if (Reads != 4) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } else { Reads = sscanf( sInput,"%d%d%lf", &Row, &Col, &Real ); if (Reads != 3) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } if(Row < 0 OR Col < 0) { fprintf(stderr, "%s: index not positive in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } if(Row > Size OR Col > Size) { if (NOT ExpansionWarningGiven) { fprintf( stderr, "%s: computed and given matrix size differ in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); ExpansionWarningGiven = YES; } Size = MAX(Row, Col); EnlargeVectors( Size, NO ); } pElement = spGetElement( Matrix, Row, Col ); if (pElement == NULL) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } pInitInfo = (ComplexNumber *)spGetInitInfo( pElement ); if (pInitInfo == NULL) { pValue[--Count].Real = Real; pValue[Count].Imag = Imag; spInstallInitInfo( pElement, (char *)(pValue + Count) ); } else { pInitInfo->Real += Real; pInitInfo->Imag += Imag; } /* Save into RHS vector if in desired column. */ if (Col == RHS_Col) { if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS RHS[Row] = Real; iRHS[Row] = Imag; #else cRHS[Row].Real = Real; cRHS[Row].Imag = Imag; #endif } else RHS[Row] = Real; } } while (Row != 0 AND Col != 0); Size = spGetSize( Matrix, YES ); if (Error = spError( Matrix ) != spOKAY) { PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) return 1; } /* Read RHS vector. */ if (NOT UseColumnAsRHS AND (NULL != fgets( sInput, BUFSIZ, pMatrixFile ))) { /* RHS vector exists, read it. */ LineNumber++; for (I = 1; I <= Size; I++) { if (I != 1 OR (strncmp( sInput, "Beginning", 8 ) == 0)) { LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } if (Complex) { #if spSEPARATED_COMPLEX_VECTORS Reads = sscanf( sInput,"%lf%lf", &RHS[I], &iRHS[I] ); #else Reads = sscanf( sInput, "%lf%lf", &cRHS[I].Real, &cRHS[I].Imag ); #endif if (Reads != 2) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } else /* Not complex. */ { Reads = sscanf( sInput, "%lf", &RHS[I] ); if (Reads != 1) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } } if (RealAsComplex AND NOT Complex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) iRHS[I] = 0.0; #else for (I = Size; I > 0; I--) { cRHS[I].Real = RHS[I]; cRHS[I].Imag = 0.0; } #endif } } /* Print out the size and the type of the matrix. */ if (NOT SolutionOnly) { IntSize = spGetSize( Matrix, NO ); printf("Matrix is %d x %d ", IntSize, IntSize); if (IntSize != Size) printf("(external size is %d x %d) ", Size, Size); if (Complex OR RealAsComplex) printf("and complex.\n",Size,Size); else printf("and real.\n",Size,Size); } if (Complex OR RealAsComplex) spSetComplex( Matrix ); else spSetReal( Matrix ); (void)fclose(pMatrixFile); return 0; } /* * INITIALIZE MATRIX ELEMENT * * This function copys the InitInfo to the Real and Imag matrix element * values. * * >>> Returned: * A zero is returns, signifying that no error can be made. */ /*ARGSUSED*/ static int Init( pElement, pInitInfo, Row, Col) RealNumber *pElement; char *pInitInfo; int Row, Col; { /* Begin `Init'. */ *pElement = *(RealNumber *)pInitInfo; if (Complex OR RealAsComplex) *(pElement+1) = *((RealNumber *)pInitInfo+1); return 0; } /* * COMPLEX ARRAY ALLOCATION * * These routines are used to check out and in arrays of complex numbers. */ static struct Bin { ComplexVector Array; struct Bin *Next; } *FirstArray, *CurrentArray = NULL; static ComplexVector CheckOutComplexArray( Count ) int Count; { struct Bin Bin, *pBin; ComplexVector Temp; /* Begin `CheckOutComplexArray'. */ if (CurrentArray == NULL OR CurrentArray->Next == NULL) { pBin = ALLOC( struct Bin, 1); Temp = ALLOC( ComplexNumber, Count ); if (pBin == NULL OR Temp == NULL) { fprintf( stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } Bin.Array = Temp; Bin.Next = NULL; *pBin = Bin; if (CurrentArray == NULL) FirstArray = CurrentArray = pBin; else { CurrentArray->Next = pBin; CurrentArray = pBin; } } else { pBin = CurrentArray; CurrentArray = pBin->Next; } return pBin->Array; } static CheckInAllComplexArrays() { /* Begin `CheckInAllComplexArrays'. */ if (CurrentArray != NULL) CurrentArray = FirstArray; return; } /* * PRINT ERROR MESSAGE * * Prints error message for matrix routines if one is required. * * >>> Argument: * Error (int) * The error status of the matrix routines. */ static PrintMatrixErrorMessage( Error ) int Error; { int Row, Col; /* Begin `PrintMatrixErrorMessage'. */ if (Error == spOKAY) return; if (Error >= spFATAL) { fprintf(stderr, "%s: fatal error detected in file `%s'.\n", ProgramName, FileName ); } else fprintf(stderr, "%s: warning in `%s'.\n", ProgramName, FileName ); if (Error == spNO_MEMORY) fprintf(stderr, " Insufficient memory available.\n"); else if (Error == spSINGULAR) { spWhereSingular( Matrix, &Row, &Col ); printf(" Singular matrix (detected at row %d, column %d).\n", Row, Col ); } else if (Error == spZERO_DIAG) { spWhereSingular( Matrix, &Row, &Col ); printf(" Zero diagonal detected at row %d, column %d.\n", Row, Col ); } else if (Error == spPANIC) fprintf(stderr, " Matrix routines being used improperly.\n"); else if (Error == spSMALL_PIVOT) fprintf(stderr, " A pivot was chosen that is smaller than the absolute threshold.\n"); return; } /* * INTERPRET THE COMMAND LINE ARGUMENTS */ static int InterpretCommandLine( ArgCount, Args ) int ArgCount; char *Args[]; { int I, FileCount = 0; char *GetProgramName(); /* Begin `InterpretCommandLine'. */ /* Determine the name of the program. */ ProgramName = GetProgramName(Args[0]); /* Step through the argument list, interpreting and deleting the options. */ for (I = 1; I < ArgCount; I++) { if (!strcmp(Args[I], "-a")) { if (ArgCount == I OR (sscanf(Args[I+1],"%lf", &AbsThreshold) != 1)) { AbsThreshold = 0.0; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-r")) { if (ArgCount == I OR (sscanf(Args[I+1],"%lf", &RelThreshold) != 1)) { RelThreshold = 0.0; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-x")) { #if spCOMPLEX RealAsComplex = YES; #else fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n", ProgramName); fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-s")) SolutionOnly = YES; else if (!strcmp(Args[I], "-c")) DiagPivoting = NO; else if (!strcmp(Args[I], "-t")) { #if TRANSPOSE Transposed = YES; #else fprintf(stderr, "%s: Sparse is not configured to solve transposed system.\n", ProgramName); fprintf(stderr," Enable TRANSPOSE in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-n")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &PrintLimit) != 1)) { PrintLimit = PRINT_LIMIT; Usage(Args[I]); } else { PrintLimitSet = YES; I++; } } else if (!strcmp(Args[I], "-i")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &Iterations) != 1)) { Iterations = 1; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-b")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &ColumnAsRHS) != 1)) { ColumnAsRHS = 1; UseColumnAsRHS = YES; } else { UseColumnAsRHS = YES; ColumnAsRHS = MAX( ColumnAsRHS, 1 ); I++; } } else if (!strcmp(Args[I], "-p")) { #if DOCUMENTATION CreatePlotFiles = YES; #else fprintf(stderr, "%s: Sparse is not configured to generate plot files.\n", ProgramName); fprintf(stderr," Enable DOCUMENTATION in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-u")) Usage( (char *)NULL ); else if (Args[I][0] == '-') Usage(Args[I]); else Args[++FileCount] = Args[I]; } return FileCount; } /* * PROGRAM NAME * * Removes path from argv[0] and returns program name. * Assumes UNIX style path names. */ static char * GetProgramName( Arg0 ) char *Arg0; { char *pTail; /* Begin `GetProgramName'. */ pTail = Arg0 + strlen(Arg0)-1; while (pTail != Arg0 AND *(pTail-1) != '/') --pTail; return pTail; } /* * USAGE WARNING * * Sends a warning to the user when the command line syntax is wrong. */ static Usage( sBadOpt ) char *sBadOpt; { /* Begin `Usage'. */ if (sBadOpt != NULL) { fprintf(stderr, "%s: unknown or deformed option `%s'.\n", ProgramName, sBadOpt); } fprintf(stderr, "\nUsage: %s [options] [file names]\n\n", ProgramName); fprintf(stderr, "Options:\n"); fprintf(stderr, " -s Print solution rather than run statistics.\n"); fprintf( stderr, " -r x Use x as relative threshold.\n"); fprintf( stderr, " -a x Use x as absolute threshold.\n"); fprintf( stderr, " -n n Print first n terms of solution vector.\n"); fprintf( stderr, " -i n Repeat build/factor/solve n times.\n"); fprintf( stderr, " -b n Use n'th column of matrix as b in Ax=b.\n"); #if DIAGONAL_PIVOTING fprintf( stderr, " -c Use complete (as opposed to diagonal) pivoting.\n"); #endif #if DOCUMENTATION fprintf( stderr, " -p Create plot files `filename.bef' and `filename.aft'.\n"); #endif #if spCOMPLEX fprintf( stderr, " -x Treat real matrix as complex with imaginary part zero.\n"); #endif #if TRANSPOSE fprintf( stderr, " -t Solve transposed system.\n"); #endif fprintf( stderr, " -u Print usage message.\n"); exit(1); } /* * ENLARGE VECTORS * * Allocate or enlarge vectors. */ static EnlargeVectors( NewSize, Clear ) int NewSize, Clear; { int I, PrevSize = MaxSize; RealVector OldRHS = RHS, iOldRHS = iRHS; ComplexVector cOldRHS = cRHS; #if spCOMPLEX # define SCALE 2 #else # define SCALE 1 #endif /* Begin `EnlargeVectors'. */ if (NewSize > PrevSize) { if (MaxSize != 0) { free( (char *)Solution ); free( (char *)RHS_Verif ); } RHS = ALLOC( RealNumber, SCALE*(NewSize+1) ); Solution = ALLOC( RealNumber, SCALE*(NewSize+1) ); RHS_Verif = ALLOC( RealNumber, SCALE*(NewSize+1) ); if (NOT RHS OR NOT Solution OR NOT RHS_Verif) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } cRHS = (ComplexVector)RHS; cSolution = (ComplexVector)Solution; cRHS_Verif = (ComplexVector)RHS_Verif; iRHS = RHS + NewSize + 1; iSolution = Solution + NewSize + 1; iRHS_Verif = RHS_Verif + NewSize + 1; /* Copy data from old RHS to new RHS. */ if (NOT Clear) { /* Copy old RHS vector to new. */ if (Complex OR RealAsComplex) { for (I = PrevSize; I > 0; I--) { #if spSEPARATED_COMPLEX_VECTORS OR LINT RHS[I] = OldRHS[I]; iRHS[I] = iOldRHS[I]; #endif #if spSEPARATED_COMPLEX_VECTORS OR LINT cRHS[I] = cOldRHS[I]; #endif } } else { for (I = PrevSize; I > 0; I--) RHS[I] = OldRHS[I]; } } if (MaxSize != 0) free( (char *)OldRHS ); MaxSize = NewSize; } /* Either completely clear or clear remaining portion of RHS vector. */ if ((NewSize > PrevSize) OR Clear) { if (Clear) { NewSize = MAX( NewSize, PrevSize ); PrevSize = 0; } if (Complex OR RealAsComplex) { for (I = NewSize; I > PrevSize; I--) { #if spSEPARATED_COMPLEX_VECTORS OR LINT RHS[I] = 0.0; iRHS[I] = 0.0; #endif #if NOT spSEPARATED_COMPLEX_VECTORS OR LINT cRHS[I].Real = 0.0; cRHS[I].Imag = 0.0; #endif } } else { for (I = NewSize; I > PrevSize; I--) RHS[I] = 0.0; } } return; }