00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 #include <stdio.h>
00049 #include <math.h>
00050 #include <ctype.h>
00051 #include <string.h>
00052 #include <time.h>
00053 #define spINSIDE_SPARSE
00054 #include "spConfig.h"
00055 #undef spINSIDE_SPARSE
00056 #include "spmatrix.h"
00057
00058 static struct Bin {
00059 ComplexVector Array;
00060 struct Bin *Next;
00061 } *FirstArray, *CurrentArray = NULL;
00062
00063
00064
00065
00066
00067 static void CheckInAllComplexArrays()
00068 {
00069 if (CurrentArray != (void*)0)
00070 CurrentArray = FirstArray;
00071 return;
00072 }
00073
00074 static int ReadMatrixFromFile();
00075 static void PrintMatrixErrorMessage(int Error);
00076 static void EnlargeVectors(int NewSize, int Clear );
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 #define PRINT_LIMIT 9
00114
00115 static char *ProgramName, *FileName;
00116 static int Size, MaxSize = 0, PrintLimit = PRINT_LIMIT, Iterations = 1, ColumnAsRHS = 1;
00117 static BOOLEAN Complex, SolutionOnly = 0, RealAsComplex = 0, Transposed = 0;
00118 static BOOLEAN CreatePlotFiles = 0, UseColumnAsRHS = 0, PrintLimitSet = 0;
00119 static BOOLEAN ExpansionWarningGiven, DiagPivoting = DIAGONAL_PIVOTING;
00120 static char Message[BUFSIZ];
00121 MatrixPtr Matrix = NULL;
00122 static double AbsThreshold = 0.0, RelThreshold = 0.0;
00123 static double *RHS, *Solution, *RHS_Verif, *iRHS, *iSolution, *iRHS_Verif;
00124 static ComplexVector cRHS, cSolution, cRHS_Verif;
00125 int Init();
00126 static void Usage(char *sBadOpt)
00127 {
00128 if (sBadOpt != NULL)
00129 fprintf(stderr, "%s: unknown or deformed option `%s'.\n", ProgramName, sBadOpt);
00130 fprintf(stderr, "\nUsage: %s [options] [file names]\n\nOptions:\n", ProgramName);
00131 fprintf(stderr, " -s Print solution rather than run statistics.\n");
00132 fprintf(stderr, " -r x Use x as relative threshold.\n");
00133 fprintf(stderr, " -a x Use x as absolute threshold.\n");
00134 fprintf(stderr, " -n n Print first n terms of solution vector.\n");
00135 fprintf(stderr, " -i n Repeat build/factor/solve n times.\n");
00136 fprintf(stderr, " -b n Use n'th column of matrix as b in Ax=b.\n");
00137 #if DIAGONAL_PIVOTING
00138 fprintf(stderr, " -c Use complete (as opposed to diagonal) pivoting.\n");
00139 #endif
00140 #if DOCUMENTATION
00141 fprintf(stderr, " -p Create plot files `filename.bef' and `filename.aft'.\n");
00142 #endif
00143 #if spCOMPLEX
00144 fprintf(stderr, " -x Treat real matrix as complex with imaginary part zero.\n");
00145 #endif
00146 #if TRANSPOSE
00147 fprintf(stderr, " -t Solve transposed system.\n");
00148 #endif
00149 fprintf(stderr, " -u Print usage message.\n");
00150 exit(1);
00151 }
00152
00153 static int InterpretCommandLine(int ArgCount, char* Args[])
00154 {
00155 int I, FileCount = 0;
00156 char *GetProgramName();
00157
00158
00159
00160
00161 char* ProgramName = GetProgramName(Args[0]);
00162
00163
00164 for (I = 1; I < ArgCount; I++)
00165 { if (!strcmp(Args[I], "-a"))
00166 {
00167 if (ArgCount == I || (sscanf(Args[I+1],"%lf", &AbsThreshold) != 1))
00168 {
00169 AbsThreshold = 0.0;
00170 Usage(Args[I]);
00171 }
00172 else I++;
00173 }
00174 else if (!strcmp(Args[I], "-r"))
00175 {
00176 if (ArgCount == I || (sscanf(Args[I+1], "%lf", &RelThreshold) != 1))
00177 {
00178 RelThreshold = 0.0;
00179 Usage(Args[I]);
00180 }
00181 else I++;
00182 }
00183 else if (!strcmp(Args[I], "-x"))
00184 {
00185 #if spCOMPLEX
00186 RealAsComplex = 1;
00187 #else
00188 fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n",
00189 ProgramName);
00190 fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n");
00191 #endif
00192 }
00193 else if (!strcmp(Args[I], "-s"))
00194 SolutionOnly = 1;
00195 else if (!strcmp(Args[I], "-c"))
00196 DiagPivoting = 0;
00197 else if (!strcmp(Args[I], "-t"))
00198 {
00199 #if TRANSPOSE
00200 Transposed = 1;
00201 #else
00202 fprintf(stderr, "%s: Sparse is not configured to solve transposed system.\n",
00203 ProgramName);
00204 fprintf(stderr," Enable TRANSPOSE in `spConfig.h'.\n");
00205 #endif
00206 }
00207 else if (!strcmp(Args[I], "-n"))
00208 {
00209 if (ArgCount == I || (sscanf(Args[I+1],"%ld", &PrintLimit) != 1))
00210 {
00211 PrintLimit = PRINT_LIMIT;
00212 Usage(Args[I]);
00213 }
00214 else
00215 {
00216 PrintLimitSet = 1; I++;
00217 }
00218 }
00219 else if (!strcmp(Args[I], "-i"))
00220 {
00221 if (ArgCount == I || (sscanf(Args[I+1],"%ld", &Iterations) != 1))
00222 {
00223 Iterations = 1; Usage(Args[I]);
00224 }
00225 else I++;
00226 }
00227 else if (!strcmp(Args[I], "-b"))
00228 {
00229 if (ArgCount == I || (sscanf(Args[I+1],"%ld", &ColumnAsRHS) != 1))
00230 {
00231 ColumnAsRHS = 1;
00232 UseColumnAsRHS = 1;
00233 }
00234 else
00235 {
00236 UseColumnAsRHS = 1;
00237 ColumnAsRHS = MAX( ColumnAsRHS, 1 );
00238 I++;
00239 }
00240 }
00241 else if (!strcmp(Args[I], "-p"))
00242 {
00243 #if DOCUMENTATION
00244 CreatePlotFiles = 1;
00245 #else
00246 fprintf(stderr, "%s: Sparse is not configured to generate plot files.\n",
00247 ProgramName);
00248 fprintf(stderr," Enable DOCUMENTATION in `spConfig.h'.\n");
00249 #endif
00250 }
00251 else if (!strcmp(Args[I], "-u"))
00252 Usage( (char *)NULL );
00253 else if (Args[I][0] == '-')
00254 Usage(Args[I]);
00255 else Args[++FileCount] = Args[I];
00256 }
00257 return FileCount;
00258 }
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312 int main(int ArgCount, char* Args[])
00313 {
00314 long I;
00315 int Error, Last, Elements, Fillins, Exponent;
00316 double MaxRHS, Residual, Determinant, iDeterminant,
00317 ConditionNumber, PseudoCondition,
00318 LargestBefore, LargestAfter, Roundoff, InfNorm,
00319 BuildTime, FactorTime, SolveTime, PartitionTime,
00320 InitialFactorTime, ConditionTime, DeterminantTime;
00321 BOOLEAN StandardInput;
00322 char j, PlotFile[BUFSIZ], ErrMsg[BUFSIZ];
00323 clock_t StartTime, BeginTime;
00324
00325 BeginTime = clock();
00326 ArgCount = InterpretCommandLine( ArgCount, Args );
00327
00328
00329
00330
00331
00332
00333 do
00334 {
00335
00336 BuildTime = FactorTime = SolveTime = 0.0;
00337 ExpansionWarningGiven = 0;
00338
00339 Matrix = spCreate(0, spCOMPLEX, &Error );
00340 if (Matrix == NULL)
00341 {
00342 fprintf(stderr, "%s: insufficient memory available.\n", ProgramName);
00343 exit(1);
00344 }
00345 if( Error >= spFATAL ) goto End;
00346
00347
00348 if (ArgCount == 0)
00349 FileName = NULL;
00350 else
00351 FileName = *(++Args);
00352 Error = ReadMatrixFromFile();
00353 if (Error) goto End;
00354 StandardInput = (FileName == NULL);
00355
00356
00357 if (spGetSize(Matrix, 1) != spGetSize(Matrix, 0))
00358 { if (Complex || RealAsComplex)
00359 { for (I = Size; I > 0; I--)
00360 cSolution[I].Real = cSolution[I].Imag = 0.0;
00361 }
00362 else
00363 { for (I = Size; I > 0; I--)
00364 Solution[I] = 0.0;
00365 }
00366 }
00367
00368
00369 (void)spInitialize(Matrix, Init);
00370
00371 #if MODIFIED_NODAL
00372 spMNA_Preorder( Matrix );
00373 #endif
00374 #if DOCUMENTATION
00375 if (CreatePlotFiles)
00376 { if (StandardInput)
00377 (void)sprintf(PlotFile, "bef");
00378 else
00379 (void)sprintf(PlotFile, "%s.bef", FileName);
00380 if (! spFileMatrix( Matrix, PlotFile, FileName, 0, 0, 0 ))
00381 { (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile);
00382 perror( ErrMsg );
00383 }
00384 }
00385 #if 0
00386 spPrint( Matrix, 0 , 0 , 0 );
00387 #endif
00388 #endif
00389 #if STABILITY
00390 if (! SolutionOnly) LargestBefore = spLargestElement(Matrix);
00391 #endif
00392 #if CONDITION
00393 if (! SolutionOnly) InfNorm = spNorm(Matrix);
00394 #endif
00395
00396 StartTime = clock();
00397 Error = spOrderAndFactor(Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting);
00398 InitialFactorTime = clock() - StartTime;
00399 PrintMatrixErrorMessage( Error );
00400 if( Error >= spFATAL )
00401 goto End;
00402
00403 #if DOCUMENTATION
00404 if (CreatePlotFiles)
00405 {
00406 if (StandardInput)
00407 (void)sprintf(PlotFile, "aft");
00408 else
00409 (void)sprintf(PlotFile, "%s.aft", FileName);
00410 if (! spFileMatrix( Matrix, PlotFile, FileName, 1, 0, 0 ))
00411 {
00412 (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile);
00413 perror( ErrMsg );
00414 }
00415 }
00416 #if 0
00417 spFileStats( Matrix, FileName, "stats" );
00418 #endif
00419 #endif
00420
00421
00422
00423
00424
00425
00426 #if TRANSPOSE
00427 if (Transposed)
00428 spSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS );
00429 else
00430 #endif
00431 spSolve(Matrix, RHS, Solution IMAG_VECTORS);
00432 if (SolutionOnly)
00433 Iterations = 0;
00434 else
00435 {
00436 #if STABILITY
00437 LargestAfter = spLargestElement(Matrix);
00438 Roundoff = spRoundoff(Matrix, LargestAfter);
00439 #endif
00440 #if CONDITION
00441 StartTime = clock();
00442 ConditionNumber = spCondition(Matrix, InfNorm, &Error);
00443 ConditionTime = clock() - StartTime;
00444 PrintMatrixErrorMessage(Error);
00445 #endif
00446 #if PSEUDOCONDITION
00447 PseudoCondition = spPseudoCondition(Matrix);
00448 #endif
00449
00450 StartTime = clock();
00451 spPartition(Matrix, spDEFAULT_PARTITION);
00452 PartitionTime = clock() - StartTime;
00453 }
00454
00455
00456 for(I = 1; I <= Iterations; I++)
00457 {
00458 StartTime = clock();
00459 (void)spInitialize(Matrix, Init);
00460 BuildTime += clock() - StartTime;
00461
00462 StartTime = clock();
00463 Error = spFactor(Matrix);
00464 FactorTime += clock() - StartTime;
00465 if( Error != spOKAY ) PrintMatrixErrorMessage( Error );
00466 if( Error >= spFATAL ) goto End;
00467 StartTime = clock();
00468
00469
00470
00471
00472
00473 #if TRANSPOSE
00474 if (Transposed)
00475 spSolveTransposed(Matrix, RHS, Solution IMAG_VECTORS );
00476 else
00477 #endif
00478 spSolve( Matrix, RHS, Solution IMAG_VECTORS );
00479 SolveTime += clock() - StartTime;
00480 }
00481
00482 if (SolutionOnly)
00483 {
00484 if (PrintLimitSet)
00485 Last = MIN( PrintLimit, Size );
00486 else
00487 Last = Size;
00488 j = ' ';
00489 }
00490 else
00491 {
00492 Last = (PrintLimit > Size) ? Size : PrintLimit;
00493 if (Last > 0) printf("Solution:\n");
00494 j = 'j';
00495 }
00496
00497 if (Complex || RealAsComplex)
00498 {
00499 #if spSEPARATED_COMPLEX_VECTORS
00500 for (I = 1; I <= Last; I++)
00501 printf("%-16.9lg %-.9lg%c\n", (double)Solution[I], (double)iSolution[I], j);
00502 #else
00503 for (I = 1; I <= Last; I++)
00504 printf("%-16.9lg %-.9lg%c\n", (double)cSolution[I].Real, (double)cSolution[I].Imag, j);
00505 #endif
00506 }
00507 else
00508 for (I = 1; I <= Last; I++)
00509 printf("%-.9lg\n", (double)Solution[I]);
00510 if (Last < Size && Last != 0)
00511 printf("Solution list truncated.\n");
00512 printf("\n");
00513
00514 #if DETERMINANT
00515
00516 if (!SolutionOnly)
00517 {
00518 StartTime = clock();
00519 #if spCOMPLEX
00520 spDeterminant( Matrix, &Exponent, &Determinant, &iDeterminant );
00521 #else
00522 spDeterminant( Matrix, &Exponent, &Determinant );
00523 #endif
00524 DeterminantTime = clock() - StartTime;
00525 if (Complex || RealAsComplex)
00526 {
00527 Determinant = hypot(Determinant, iDeterminant);
00528 while (Determinant >= 10.0)
00529 {
00530 Determinant *= 0.1; Exponent++;
00531 }
00532 }
00533 }
00534 #else
00535 Determinant = 0.0; Exponent = 0;
00536 #endif
00537
00538 #if MULTIPLICATION
00539 if (!SolutionOnly)
00540 {
00541
00542
00543 MaxRHS = 0.0;
00544 if (Complex || RealAsComplex)
00545 {
00546 #if spSEPARATED_COMPLEX_VECTORS
00547 for (I = 1; I <= Size; I++)
00548 {
00549 if (ABS(RHS[I]) > MaxRHS)
00550 MaxRHS = ABS(RHS[I]);
00551 if (ABS(iRHS[I]) > MaxRHS)
00552 MaxRHS = ABS(iRHS[I]);
00553 }
00554 #else
00555 for (I = 1; I <= Size; I++)
00556 {
00557 if (ABS(cRHS[I].Real) > MaxRHS)
00558 MaxRHS = ABS(cRHS[I].Real);
00559 if (ABS(cRHS[I].Imag) > MaxRHS)
00560 MaxRHS = ABS(cRHS[I].Imag);
00561 }
00562 #endif
00563 }
00564 else
00565 for (I = 1; I <= Size; I++)
00566 if (ABS(RHS[I]) > MaxRHS)
00567 MaxRHS = ABS(RHS[I]);
00568
00569
00570 (void)spInitialize(Matrix, Init);
00571 #if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS
00572 if (Transposed)
00573 spMultTransposed(Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution);
00574 else
00575 spMultiply(Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution);
00576 #else
00577 if (Transposed)
00578 spMultTransposed(Matrix, RHS_Verif, Solution);
00579 else
00580 spMultiply(Matrix, RHS_Verif, Solution);
00581 #endif
00582
00583
00584 Residual = 0.0;
00585 if (Complex || RealAsComplex)
00586 {
00587 #if spSEPARATED_COMPLEX_VECTORS
00588 for (I = 1; I <= Size; I++)
00589 Residual += ABS(RHS[I] - RHS_Verif[I]) + ABS(iRHS[I] - iRHS_Verif[I]);
00590 #else
00591 for (I = 1; I <= Size; I++)
00592 Residual += ABS(cRHS[I].Real - cRHS_Verif[I].Real) + ABS(cRHS[I].Imag - cRHS_Verif[I].Imag);
00593 #endif
00594 }
00595 else
00596 for (I = 1; I <= Size; I++)
00597 Residual += ABS(RHS[I] - RHS_Verif[I]);
00598 }
00599 #endif
00600
00601
00602 if (! SolutionOnly)
00603 {
00604 Elements = spElementCount(Matrix);
00605 Fillins = spFillinCount(Matrix);
00606
00607 printf("Initial factor time = %.2lf.\n", InitialFactorTime);
00608 printf("Partition time = %.2lf.\n", PartitionTime);
00609 if (Iterations > 0)
00610 {
00611 printf("Build time = %.3lf.\n", (BuildTime/Iterations));
00612 printf("Factor time = %.3lf.\n",(FactorTime/Iterations));
00613 printf("Solve time = %.3lf.\n", (SolveTime/Iterations));
00614 }
00615 #if STABILITY
00616 printf("Condition time = %.2lf.\n", ConditionTime);
00617 #endif
00618 #if DETERMINANT
00619 printf("Determinant time = %.2lf.\n", DeterminantTime);
00620 #endif
00621 printf("\nTotal number of elements = %d.\n", Elements);
00622 printf("Average number of elements per row initially = %.2lf.\n",
00623 (double)(Elements - Fillins) / (double)spGetSize(Matrix, 0));
00624 printf("Total number of fill-ins = %d.\n", Fillins);
00625 #if DETERMINANT || MULTIPLICATION || PSEUDOCONDITION || CONDITION || STABILITY
00626 putchar('\n');
00627 #endif
00628 #if STABILITY
00629 if (LargestBefore != 0.0)
00630 printf("Growth = %.2lg.\n", LargestAfter / LargestBefore);
00631 printf("Max error in matrix = %.2lg.\n", Roundoff);
00632 #endif
00633 #if STABILITY
00634 if(ABS(ConditionNumber) > 10 * DBL_MIN);
00635 printf("Condition number = %.2lg.\n", 1.0 / ConditionNumber);
00636 #endif
00637 #if CONDITION && STABILITY
00638 printf("Estimated upper bound of error in solution = %.2lg.\n",
00639 Roundoff / ConditionNumber);
00640 #endif
00641 #if PSEUDOCONDITION
00642 printf("PseudoCondition = %.2lg.\n", PseudoCondition);
00643 #endif
00644 #if DETERMINANT
00645 printf("Determinant = %.3lg", (double)Determinant );
00646 if (Determinant != 0.0 && Exponent != 0)
00647 printf("e%d", Exponent);
00648 putchar('.'); putchar('\n');
00649 #endif
00650 #if MULTIPLICATION
00651 if (MaxRHS != 0.0)
00652 printf("Normalized residual = %.2lg.\n", (Residual / MaxRHS));
00653 #endif
00654 }
00655
00656 End:;
00657 (void)fflush(stdout);
00658 CheckInAllComplexArrays();
00659 spDestroy(Matrix);
00660 Matrix = NULL;
00661 } while( --ArgCount > 0 );
00662
00663 if (!SolutionOnly)
00664 {
00665 printf("\nAggregate resource usage:\n");
00666 printf(" Time required = %.2lf seconds.\n", (double)(clock() - BeginTime) / CLOCKS_PER_SEC );
00667
00668 }
00669 return 0;
00670 }
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 static int ReadMatrixFromFile()
00700 {
00701 long I, Reads;
00702 FILE *pMatrixFile;
00703 char sInput[BUFSIZ], sType[BUFSIZ], *p;
00704 int Error, Row, Col, Count = 0, LineNumber, RHS_Col, IntSize;
00705 double Real, Imag = 0.0;
00706 double* pElement;
00707 ComplexNumber *pValue, *pInitInfo, *CheckOutComplexArray();
00708 static char *EndOfFile = "%s: unexpected end of file `%s' at line %d.\n",
00709 *Syntax = "%s: syntax error in file `%s' at line %d.\n";
00710
00711
00712 if (!SolutionOnly) putchar('\n');
00713
00714 if (FileName == NULL)
00715 {
00716 FileName = "standard input";
00717 pMatrixFile = stdin;
00718 }
00719 else
00720 {
00721 pMatrixFile = fopen(FileName, "r");
00722 if (pMatrixFile == NULL)
00723 {
00724 fprintf(stderr, "%s: file %s was not found.\n", ProgramName, FileName);
00725 return 1;
00726 }
00727 }
00728 Complex = 0;
00729 LineNumber = 1;
00730
00731
00732 if (NULL == fgets(Message, BUFSIZ, pMatrixFile))
00733 {
00734 fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1;
00735 }
00736
00737
00738 if (!strncmp( Message, "Starting", 8 ))
00739 {
00740 if (strncmp( Message, "Starting complex", 15 ) == 0)
00741 Complex = 1;
00742 LineNumber++;
00743 if (NULL == fgets( Message, BUFSIZ, pMatrixFile ))
00744 {
00745 fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber);
00746 return 1;
00747 }
00748 }
00749 if (!SolutionOnly) printf("%-s\n", Message);
00750
00751
00752 LineNumber++;
00753 if (NULL == fgets( sInput, BUFSIZ, pMatrixFile ))
00754 {
00755 fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1;
00756 }
00757 if ((Reads = sscanf( sInput,"%d %s", &Size, sType )) < 1)
00758 {
00759 fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1;
00760 }
00761 if (Reads == 2)
00762 {
00763 for (p = sType; *p != '\0'; p++)
00764 if (isupper(*p)) *p += 'a'-'A';
00765 if (strncmp( sType, "complex", 7 ) == 0)
00766 Complex = 1;
00767 else if (strncmp( sType, "real", 7 ) == 0)
00768 Complex = 0;
00769 else
00770 {
00771 fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1;
00772 }
00773 }
00774 EnlargeVectors( Size, 1 );
00775 RHS_Col = MIN( Size, ColumnAsRHS );
00776
00777 #if ! spCOMPLEX
00778 if (Complex)
00779 {
00780 fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n",
00781 ProgramName);
00782 fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n");
00783 return 1;
00784 }
00785 #endif
00786 #if ! REAL
00787 if (!(Complex || RealAsComplex))
00788 {
00789 fprintf(stderr, "%s: Sparse is not configured to solve real matrices.\n",
00790 ProgramName);
00791 fprintf(stderr," Enable REAL in `spConfig.h'.\n");
00792 return 1;
00793 }
00794 #endif
00795
00796
00797 do
00798 { if (Count == 0)
00799 pValue = CheckOutComplexArray(Count = 1000);
00800 LineNumber++;
00801 if (NULL == fgets( sInput, BUFSIZ, pMatrixFile ))
00802 {
00803 fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n",
00804 ProgramName, FileName, LineNumber);
00805 return 1;
00806 }
00807 if (Complex)
00808 {
00809 Reads = sscanf( sInput,"%d%d%lf%lf", &Row, &Col, &Real, &Imag );
00810 if (Reads != 4)
00811 { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n",
00812 ProgramName, FileName, LineNumber);
00813 return 1;
00814 }
00815 }
00816 else
00817 {
00818 Reads = sscanf( sInput,"%d%d%lf", &Row, &Col, &Real );
00819 if (Reads != 3)
00820 {
00821 fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n",
00822 ProgramName, FileName, LineNumber);
00823 return 1;
00824 }
00825 }
00826 if(Row < 0 || Col < 0)
00827 {
00828 fprintf(stderr, "%s: index not positive in file `%s' at line %d.\n",
00829 ProgramName, FileName, LineNumber);
00830 return 1;
00831 }
00832 if(Row > Size || Col > Size)
00833 {
00834 if (!ExpansionWarningGiven)
00835 {
00836 fprintf( stderr,
00837 "%s: computed and given matrix size differ in file `%s' at line %d.\n",
00838 ProgramName, FileName, LineNumber);
00839 ExpansionWarningGiven = 1;
00840 }
00841 Size = MAX(Row, Col);
00842 EnlargeVectors( Size, 0 );
00843 }
00844 pElement = spGetElement(Matrix, Row, Col);
00845 if (pElement == NULL)
00846 {
00847 fprintf(stderr, "%s: insufficient memory available.\n", ProgramName);
00848 exit(1);
00849 }
00850 pInitInfo = (ComplexNumber*)spGetInitInfo(pElement);
00851 if (pInitInfo == NULL)
00852 {
00853 pValue[--Count].Real = Real;
00854 pValue[Count].Imag = Imag;
00855 spInstallInitInfo(pElement, (char *)(pValue + Count));
00856 }
00857 else
00858 {
00859 pInitInfo->Real += Real;
00860 pInitInfo->Imag += Imag;
00861 }
00862
00863
00864 if (Col == RHS_Col)
00865 {
00866 if (Complex || RealAsComplex)
00867 {
00868 #if spSEPARATED_COMPLEX_VECTORS
00869 RHS[Row] = Real;
00870 iRHS[Row] = Imag;
00871 #else
00872 cRHS[Row].Real = Real;
00873 cRHS[Row].Imag = Imag;
00874 #endif
00875 }
00876 else RHS[Row] = Real;
00877 }
00878 } while (Row != 0 && Col != 0);
00879
00880 Size = spGetSize( Matrix, 1 );
00881 if (Error = spError( Matrix ) != spOKAY)
00882 {
00883 PrintMatrixErrorMessage( Error );
00884 if(Error >= spFATAL) return 1;
00885 }
00886
00887
00888 if (! UseColumnAsRHS && (NULL != fgets( sInput, BUFSIZ, pMatrixFile )))
00889 {
00890
00891 LineNumber++;
00892 for (I = 1; I <= Size; I++)
00893 {
00894 if (I != 1 || (strncmp( sInput, "Beginning", 8 ) == 0))
00895 {
00896 LineNumber++;
00897 if (NULL == fgets( sInput, BUFSIZ, pMatrixFile ))
00898 {
00899 fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n",
00900 ProgramName, FileName, LineNumber);
00901 return 1;
00902 }
00903 }
00904 if (Complex)
00905 {
00906 #if spSEPARATED_COMPLEX_VECTORS
00907 Reads = sscanf( sInput,"%lf%lf", &RHS[I], &iRHS[I] );
00908 #else
00909 Reads = sscanf( sInput, "%lf%lf", &cRHS[I].Real, &cRHS[I].Imag );
00910 #endif
00911 if (Reads != 2)
00912 { fprintf(stderr,
00913 "%s: syntax error in file `%s' at line %d.\n",
00914 ProgramName, FileName, LineNumber);
00915 return 1;
00916 }
00917 }
00918 else
00919 {
00920 Reads = sscanf( sInput, "%lf", &RHS[I] );
00921 if (Reads != 1)
00922 {
00923 fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n",
00924 ProgramName, FileName, LineNumber);
00925 return 1;
00926 }
00927 }
00928 }
00929 if (RealAsComplex && ! Complex)
00930 {
00931 #if spSEPARATED_COMPLEX_VECTORS
00932 for (I = 1; I <= Size; I++) iRHS[I] = 0.0;
00933 #else
00934 for (I = Size; I > 0; I--)
00935 {
00936 cRHS[I].Real = RHS[I];
00937 cRHS[I].Imag = 0.0;
00938 }
00939 #endif
00940 }
00941 }
00942
00943
00944 if (! SolutionOnly)
00945 {
00946 IntSize = spGetSize( Matrix, 0 );
00947 printf("Matrix is %d x %d ", IntSize, IntSize);
00948 if (IntSize != Size)
00949 printf("(external size is %d x %d) ", Size, Size);
00950 if (Complex || RealAsComplex)
00951 printf("and complex.\n",Size,Size);
00952 else
00953 printf("and real.\n",Size,Size);
00954 }
00955 if (Complex || RealAsComplex)
00956 spSetComplex(Matrix);
00957 else
00958 spSetReal(Matrix);
00959 (void)fclose(pMatrixFile);
00960 return 0;
00961 }
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974 static int Init(double *pElement, char *pInitInfo)
00975 {
00976 *pElement = *(double*)pInitInfo;
00977 if (Complex || RealAsComplex)
00978 *(pElement+1) = *(double *)pInitInfo+1;
00979 return 0;
00980 }
00981
00982 static ComplexVector CheckOutComplexArray(int Count)
00983 {
00984 struct Bin Bin, *pBin;
00985 ComplexVector Temp;
00986
00987 if (CurrentArray == NULL || CurrentArray->Next == NULL)
00988 {
00989 pBin = ALLOC( struct Bin, 1);
00990 Temp = ALLOC( ComplexNumber, Count );
00991 if (pBin == NULL || Temp == NULL)
00992 {
00993 fprintf( stderr, "%s: insufficient memory available.\n", ProgramName);
00994 exit(1);
00995 }
00996 Bin.Array = Temp; Bin.Next = NULL;
00997 *pBin = Bin;
00998 if (CurrentArray == NULL)
00999 FirstArray = CurrentArray = pBin;
01000 else
01001 CurrentArray->Next = CurrentArray = pBin;
01002 }
01003 else
01004 {
01005 pBin = CurrentArray;
01006 CurrentArray = pBin->Next;
01007 }
01008 return pBin->Array;
01009 }
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020 static void PrintMatrixErrorMessage(int Error)
01021 {
01022 int Row, Col;
01023
01024 if (Error == spOKAY)
01025 return;
01026 if (Error >= spFATAL)
01027 fprintf(stderr, "%s: fatal error detected in file `%s'.\n", ProgramName, FileName );
01028 else
01029 fprintf(stderr, "%s: warning in `%s'.\n", ProgramName, FileName );
01030 if (Error == spNO_MEMORY)
01031 fprintf(stderr, " Insufficient memory available.\n");
01032 else if (Error == spSINGULAR)
01033 {
01034 spWhereSingular( Matrix, &Row, &Col );
01035 printf(" Singular matrix (detected at row %d, column %d).\n", Row, Col );
01036 }
01037 else if (Error == spZERO_DIAG)
01038 {
01039 spWhereSingular( Matrix, &Row, &Col );
01040 printf(" Zero diagonal detected at row %d, column %d.\n", Row, Col );
01041 }
01042 else if (Error == spPANIC)
01043 fprintf(stderr, " Matrix routines being used improperly.\n");
01044 else if (Error == spSMALL_PIVOT)
01045 fprintf(stderr, " A pivot was chosen that is smaller than the absolute threshold.\n");
01046 return;
01047 }
01048
01049
01050
01051
01052
01053
01054
01055 static char* GetProgramName(char *Arg0)
01056 {
01057 char *pTail = Arg0 + strlen(Arg0)-1;
01058 while (pTail != Arg0 && *(pTail-1) != '/')
01059 --pTail;
01060 return pTail;
01061 }
01062
01063
01064
01065
01066
01067
01068 static void EnlargeVectors(int NewSize, int Clear)
01069 {
01070 int I, PrevSize = MaxSize;
01071 double* OldRHS = RHS, *iOldRHS = iRHS;
01072 ComplexVector cOldRHS = cRHS;
01073
01074
01075
01076
01077
01078 #define SCALE 1
01079
01080 if (NewSize > PrevSize)
01081 {
01082 if (MaxSize != 0)
01083 {
01084 free( (char *)Solution );
01085 free( (char *)RHS_Verif );
01086 }
01087 RHS = ALLOC( double, SCALE*(NewSize+1) );
01088 Solution = ALLOC( double, SCALE*(NewSize+1) );
01089 RHS_Verif = ALLOC( double, SCALE*(NewSize+1) );
01090 if (! RHS || ! Solution || ! RHS_Verif)
01091 {
01092 fprintf(stderr, "%s: insufficient memory available.\n", ProgramName);
01093 exit(1);
01094 }
01095 cRHS = (ComplexVector)RHS;
01096 cSolution = (ComplexVector)Solution;
01097 cRHS_Verif = (ComplexVector)RHS_Verif;
01098 iRHS = RHS + NewSize + 1;
01099 iSolution = Solution + NewSize + 1;
01100 iRHS_Verif = RHS_Verif + NewSize + 1;
01101
01102
01103 if (!Clear)
01104
01105 if (Complex || RealAsComplex)
01106 {
01107 for (I = PrevSize; I > 0; I--)
01108 {
01109 #if spSEPARATED_COMPLEX_VECTORS || LINT
01110 RHS[I] = OldRHS[I];
01111 iRHS[I] = iOldRHS[I];
01112 #endif
01113 #if spSEPARATED_COMPLEX_VECTORS || LINT
01114 cRHS[I] = cOldRHS[I];
01115 #endif
01116 }
01117 }
01118 else
01119 for (I = PrevSize; I > 0; I--)
01120 RHS[I] = OldRHS[I];
01121 if (MaxSize != 0) free( (char *)OldRHS );
01122 MaxSize = NewSize;
01123 }
01124
01125
01126 if ((NewSize > PrevSize) || Clear)
01127 { if (Clear)
01128 {
01129 NewSize = MAX( NewSize, PrevSize );
01130 PrevSize = 0;
01131 }
01132 if (Complex || RealAsComplex)
01133 for (I = NewSize; I > PrevSize; I--)
01134 {
01135 #if spSEPARATED_COMPLEX_VECTORS || LINT
01136 RHS[I] = 0.0;
01137 iRHS[I] = 0.0;
01138 #endif
01139 #if ! spSEPARATED_COMPLEX_VECTORS || LINT
01140 cRHS[I].Real = 0.0;
01141 cRHS[I].Imag = 0.0;
01142 #endif
01143 }
01144 else
01145 for (I = NewSize; I > PrevSize; I--)
01146 RHS[I] = 0.0;
01147 }
01148 return;
01149 }
01150