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 #define spINSIDE_SPARSE
00033 #include "spconfig.h"
00034 #include "spmatrix.h"
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
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 void spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
00105 {
00106 int J = 0, I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
00107 int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
00108 double Magnitude, SmallestDiag, SmallestElement,
00109 LargestElement = 0.0, LargestDiag = 0.0;
00110 ElementPtr pElement, pImagElements[PRINTER_WIDTH/10+1];
00111
00112 ASSERT( IS_SPARSE( Matrix ) );
00113 Size = Matrix->Size;
00114
00115
00116 # if TRANSLATE
00117 Top = Matrix->AllocatedExtSize;
00118 #else
00119 Top = Matrix->AllocatedSize;
00120 #endif
00121 CALLOC( PrintOrdToIntRowMap, int, Top + 1 );
00122 CALLOC( PrintOrdToIntColMap, int, Top + 1 );
00123 if ( PrintOrdToIntRowMap == NULL || PrintOrdToIntColMap == NULL)
00124 {
00125 Matrix->Error = spNO_MEMORY; return;
00126 }
00127 for (I = 1; I <= Size; I++)
00128 {
00129 PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I;
00130 PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I;
00131 }
00132
00133
00134 for (J = 1, I = 1; I <= Top; I++)
00135 if (PrintOrdToIntRowMap[I] != 0)
00136 PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ];
00137 for (J = 1, I = 1; I <= Top; I++)
00138 if (PrintOrdToIntColMap[I] != 0)
00139 PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ];
00140
00141 if (Header)
00142 {
00143 printf("MATRIX SUMMARY\n\n");
00144 printf("Size of matrix = %1u x %1u.\n", Size, Size);
00145 if ( Matrix->Reordered && PrintReordered )
00146 printf("Matrix has been reordered.\n");
00147 putchar('\n');
00148 if ( Matrix->Factored )
00149 printf("Matrix after factorization:\n");
00150 else
00151 printf("Matrix before factorization:\n");
00152 SmallestElement = DBL_MAX;
00153 SmallestDiag = SmallestElement;
00154 }
00155
00156 Columns = PRINTER_WIDTH;
00157 if (Header) Columns -= 5;
00158 if (Data) Columns = (Columns+1) / 10;
00159
00160
00161
00162
00163
00164 J = 0;
00165 while ( J <= Size )
00166
00167
00168 {
00169 StopCol = StartCol + Columns - 1;
00170 if (StopCol > Size)
00171 StopCol = Size;
00172
00173
00174 if (Header)
00175 {
00176 if (Data)
00177 {
00178 printf(" ");
00179 for (I = StartCol; I <= StopCol; I++)
00180 {
00181 if (PrintReordered)
00182 Col = I;
00183 else
00184 Col = PrintOrdToIntColMap[I];
00185 printf(" %9d", Matrix->IntToExtColMap[ Col ]);
00186 }
00187 printf("\n\n");
00188 }
00189 else
00190 {
00191 if (PrintReordered)
00192 printf("Columns %1d to %1d.\n",StartCol,StopCol);
00193 else
00194 {
00195 printf("Columns %1d to %1d.\n",
00196 Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ],
00197 Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]);
00198 }
00199 }
00200 }
00201
00202
00203 for (I = 1; I <= Size; I++)
00204 {
00205 if (PrintReordered)
00206 Row = I;
00207 else
00208 Row = PrintOrdToIntRowMap[I];
00209
00210 if (Header)
00211 {
00212 if (PrintReordered && !Data)
00213 printf("%4d", I);
00214 else
00215 printf("%4d", Matrix->IntToExtRowMap[ Row ]);
00216 if (!Data) putchar(' ');
00217 }
00218
00219
00220 for (J = StartCol; J <= StopCol; J++)
00221 {
00222 if (PrintReordered)
00223 Col = J;
00224 else
00225 Col = PrintOrdToIntColMap[J];
00226
00227 pElement = Matrix->FirstInCol[Col];
00228 while(pElement != NULL && pElement->Row != Row)
00229 pElement = pElement->NextInCol;
00230
00231 if (Data)
00232 pImagElements[J - StartCol] = pElement;
00233
00234 if (pElement != NULL)
00235
00236
00237 {
00238 if (Data)
00239 printf(" %9.3g", (double)pElement->Real);
00240 else
00241 putchar('x');
00242
00243
00244 if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElement)
00245 LargestElement = Magnitude;
00246 if ((Magnitude < SmallestElement) && (Magnitude != 0.0))
00247 SmallestElement = Magnitude;
00248 ElementCount++;
00249 }
00250
00251
00252 else
00253 {
00254 if (Data)
00255 printf(" ...");
00256 else
00257 putchar('.');
00258 }
00259 }
00260 putchar('\n');
00261 }
00262
00263
00264 StartCol = StopCol;
00265 StartCol++;
00266 putchar('\n');
00267 }
00268 if (Header)
00269 {
00270 printf("\nLargest element in matrix = %-1.4g.\n", LargestElement);
00271 printf("Smallest element in matrix = %-1.4g.\n", SmallestElement);
00272
00273
00274 for (I = 1; I <= Size; I++)
00275 {
00276 if (Matrix->Diag[I] != NULL)
00277 {
00278 Magnitude = ELEMENT_MAG( Matrix->Diag[I] );
00279 if ( Magnitude > LargestDiag ) LargestDiag = Magnitude;
00280 if ( Magnitude < SmallestDiag ) SmallestDiag = Magnitude;
00281 }
00282 }
00283
00284
00285 if ( Matrix->Factored )
00286 {
00287 printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag);
00288 printf("Smallest diagonal element = %-1.4g.\n", SmallestDiag);
00289 }
00290 else
00291 {
00292 printf("\nLargest pivot element = %-1.4g.\n", LargestDiag);
00293 printf("Smallest pivot element = %-1.4g.\n", SmallestDiag);
00294 }
00295
00296
00297 printf("\nDensity = %2.2f%%.\n", ((double)(ElementCount * 100)) /
00298 ((double)(Size * Size)));
00299 if (!Matrix->NeedsOrdering)
00300 printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
00301 }
00302 putchar('\n');
00303 (void)fflush(stdout);
00304
00305 FREE(PrintOrdToIntColMap);
00306 FREE(PrintOrdToIntRowMap);
00307 return;
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350 int spFileMatrix(MatrixPtr Matrix, char *File, char *Label, int Reordered, int Data, int Header)
00351 {
00352 int I, Size, Row, Col, Err;
00353 ElementPtr pElement;
00354 FILE *pMatrixFile;
00355 if ((pMatrixFile = fopen(File, "w")) == NULL)
00356 return 0;
00357 Size = Matrix->Size;
00358 if (Header)
00359 {
00360 if (Matrix->Factored && Data)
00361 Err = fprintf(pMatrixFile, "Warning : The following matrix is factored in to LU form.\n");
00362 if (Err < 0) return 0;
00363 if (fprintf(pMatrixFile, "%s\n", Label) < 0) return 0;
00364 Err = fprintf( pMatrixFile, "%d\t%s\n", Size, (Matrix->Complex ? "complex" : "real"));
00365 if (Err < 0) return 0;
00366 }
00367 if (!Data)
00368 { for (I = 1; I <= Size; I++)
00369 { pElement = Matrix->FirstInCol[I];
00370 while (pElement != NULL)
00371 { if (Reordered)
00372 { Row = pElement->Row;
00373 Col = I;
00374 }
00375 else
00376 { Row = Matrix->IntToExtRowMap[pElement->Row];
00377 Col = Matrix->IntToExtColMap[I];
00378 }
00379 pElement = pElement->NextInCol;
00380 if (fprintf(pMatrixFile, "%d\t%d\n", Row, Col) < 0) return 0;
00381 }
00382 }
00383
00384 if (Header)
00385 if (fprintf(pMatrixFile, "0\t0\n") < 0) return 0;
00386 }
00387 if (Data && !Matrix->Complex)
00388 { for (I = 1; I <= Size; I++)
00389 { pElement = Matrix->FirstInCol[I];
00390 while (pElement != NULL)
00391 { Row = Matrix->IntToExtRowMap[pElement->Row];
00392 Col = Matrix->IntToExtColMap[I];
00393 Err = fprintf
00394 ( pMatrixFile,"%d\t%d\t%-.15g\n",
00395 Row, Col, (double)pElement->Real
00396 );
00397 if (Err < 0) return 0;
00398 pElement = pElement->NextInCol;
00399 }
00400 }
00401
00402 if (Header)
00403 if (fprintf(pMatrixFile,"0\t0\t0.0\n") < 0) return 0;
00404 }
00405 if (fclose(pMatrixFile) < 0) return 0;
00406 return 1;
00407 }
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 int spFileVector(MatrixPtr Matrix,char * File, double* RHS)
00446 {
00447 int I, Size;
00448 FILE *pMatrixFile;
00449
00450 ASSERT( IS_SPARSE( Matrix ) && RHS != NULL)
00451
00452
00453 if ((pMatrixFile = fopen(File,"a")) == NULL)
00454 return 0;
00455
00456
00457 #if !ARRAY_OFFSET
00458 #if spCOMPLEX
00459 if (Matrix->Complex)
00460 {
00461 #if spSEPARATED_COMPLEX_VECTORS
00462 ASSERT(iRHS != NULL)
00463 --RHS; --iRHS;
00464 #else
00465 RHS -= 2;
00466 #endif
00467 }
00468 else
00469 #endif
00470 --RHS;
00471 #endif
00472
00473
00474
00475 Size = Matrix->Size;
00476 #if REAL && spCOMPLEX
00477 else
00478 #endif
00479 #if REAL
00480 for (I = 1; I <= Size; I++)
00481 if (fprintf(pMatrixFile, "%-.15g\n", (double)RHS[I]) < 0)
00482 return 0;
00483 #endif
00484
00485
00486 if (fclose(pMatrixFile) < 0) return 0;
00487 return 1;
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 int spFileStats(MatrixPtr Matrix, char *File, char *Label)
00526 {
00527 int Size, I, NumberOfElements;
00528 ElementPtr pElement;
00529 double Data, LargestElement, SmallestElement;
00530 FILE *pStatsFile;
00531 ASSERT( IS_SPARSE( Matrix ) );
00532
00533 if ((pStatsFile = fopen(File, "a")) == NULL)
00534 return 0;
00535
00536
00537 Size = Matrix->Size;
00538 if (!Matrix->Factored)
00539 fprintf(pStatsFile, "Matrix has not been factored.\n");
00540 fprintf(pStatsFile, "||| Starting new matrix |||\n");
00541 fprintf(pStatsFile, "%s\n", Label);
00542 if (Matrix->Complex)
00543 fprintf(pStatsFile, "Matrix is complex.\n");
00544 else
00545 fprintf(pStatsFile, "Matrix is real.\n");
00546 fprintf(pStatsFile," Size = %d\n",Size);
00547
00548
00549 NumberOfElements = 0;
00550 LargestElement = 0.0;
00551 SmallestElement = DBL_MAX;
00552
00553 for (I = 1; I <= Size; I++)
00554 {
00555 pElement = Matrix->FirstInCol[I];
00556 while (pElement != NULL)
00557 {
00558 NumberOfElements++;
00559 Data = ELEMENT_MAG(pElement);
00560 if (Data > LargestElement)
00561 LargestElement = Data;
00562 if (Data < SmallestElement && Data != 0.0)
00563 SmallestElement = Data;
00564 pElement = pElement->NextInCol;
00565 }
00566 }
00567 SmallestElement = MIN( SmallestElement, LargestElement );
00568
00569 fprintf(pStatsFile, " Initial number of elements = %d\n",
00570 NumberOfElements - Matrix->Fillins);
00571 fprintf(pStatsFile,
00572 " Initial average number of elements per row = %f\n",
00573 (double)(NumberOfElements - Matrix->Fillins) / (double)Size);
00574 fprintf(pStatsFile, " Fill-ins = %d\n",Matrix->Fillins);
00575 fprintf(pStatsFile, " Average number of fill-ins per row = %f%%\n",
00576 (double)Matrix->Fillins / (double)Size);
00577 fprintf(pStatsFile, " Total number of elements = %d\n",
00578 NumberOfElements);
00579 fprintf(pStatsFile, " Average number of elements per row = %f\n",
00580 (double)NumberOfElements / (double)Size);
00581 fprintf(pStatsFile," Density = %f%%\n",
00582 (double)(100.0*NumberOfElements)/(double)(Size*Size));
00583 fprintf(pStatsFile," Relative Threshold = %e\n", Matrix->RelThreshold);
00584 fprintf(pStatsFile," Absolute Threshold = %e\n", Matrix->AbsThreshold);
00585 fprintf(pStatsFile," Largest Element = %e\n", LargestElement);
00586 fprintf(pStatsFile," Smallest Element = %e\n\n\n", SmallestElement);
00587
00588 (void)fclose(pStatsFile);
00589 return 1;
00590 }