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
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 #define spINSIDE_SPARSE
00065 #include "spconfig.h"
00066 #include "spmatrix.h"
00067
00068 #if MODIFIED_NODAL
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
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 static int CountTwins(MatrixPtr Matrix, int Col, ElementPtr *ppTwin1,ElementPtr * ppTwin2 )
00166 {
00167 int Row, Twins = 0;
00168 ElementPtr pTwin1, pTwin2;
00169
00170
00171
00172 pTwin1 = Matrix->FirstInCol[Col];
00173 while (pTwin1 != NULL)
00174 { if (ABSS(pTwin1->Real) == 1.0)
00175 { Row = pTwin1->Row;
00176 pTwin2 = Matrix->FirstInCol[Row];
00177 while ((pTwin2 != NULL) && (pTwin2->Row != Col))
00178 pTwin2 = pTwin2->NextInCol;
00179 if ((pTwin2 != NULL) && (ABSS(pTwin2->Real) == 1.0))
00180 {
00181 if (++Twins >= 2) return Twins;
00182 (*ppTwin1 = pTwin1)->Col = Col;
00183 (*ppTwin2 = pTwin2)->Col = Row;
00184 }
00185 }
00186 pTwin1 = pTwin1->NextInCol;
00187 }
00188 return Twins;
00189 }
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 static void SwapCols(MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2)
00202 {
00203 int Col1 = pTwin1->Col, Col2 = pTwin2->Col;
00204
00205 SWAP (ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
00206 SWAP (int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
00207 #if TRANSLATE
00208 Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]] = Col2;
00209 Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]] = Col1;
00210 #endif
00211
00212 Matrix->Diag[Col1] = pTwin2;
00213 Matrix->Diag[Col2] = pTwin1;
00214 Matrix->NumberOfInterchangesIsOdd = !Matrix->NumberOfInterchangesIsOdd;
00215 return;
00216 }
00217
00218 void spMNA_Preorder(MatrixPtr Matrix)
00219 {
00220 int J, Size, Twins, StartAt = 1;
00221 ElementPtr pTwin1, pTwin2;
00222 BOOLEAN Swapped, AnotherPassNeeded;
00223
00224 ASSERT( IS_VALID(Matrix) && !Matrix->Factored );
00225 if (Matrix->RowsLinked) return;
00226 Size = Matrix->Size;
00227 Matrix->Reordered = 1;
00228 do
00229 {
00230 AnotherPassNeeded = Swapped = 0;
00231
00232
00233 for (J = StartAt; J <= Size; J++)
00234 if (Matrix->Diag[J] == NULL)
00235 {
00236 Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
00237 if (Twins == 1)
00238 {
00239 SwapCols( Matrix, pTwin1, pTwin2 );
00240 Swapped = 1;
00241 }
00242 else if ((Twins > 1) && !AnotherPassNeeded)
00243 { AnotherPassNeeded = 1;
00244 StartAt = J;
00245 }
00246 }
00247
00248 if (AnotherPassNeeded)
00249 for (J = StartAt; !Swapped && (J <= Size); J++)
00250 if (Matrix->Diag[J] == NULL)
00251 {
00252 Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
00253 SwapCols( Matrix, pTwin1, pTwin2 );
00254 Swapped = 0;
00255 }
00256 } while (AnotherPassNeeded);
00257 return;
00258 }
00259 #endif
00260
00261 #if SCALING
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
00313
00314
00315
00316
00317 extern void spcLinkRows(MatrixPtr Matrix);
00318
00319
00320
00321
00322 void spScale(MatrixPtr Matrix, double * RHS_ScaleFactors, double* SolutionScaleFactors )
00323 {
00324 ElementPtr pElement;
00325 int I, lSize, *pExtOrder;
00326 double ScaleFactor;
00327
00328
00329
00330 ASSERT( IS_VALID(Matrix) && !Matrix->Factored );
00331 if (!Matrix->RowsLinked) spcLinkRows( Matrix );
00332 #if REAL
00333 lSize = Matrix->Size;
00334
00335
00336 #if !ARRAY_OFFSET
00337 --RHS_ScaleFactors;
00338 --SolutionScaleFactors;
00339 #endif
00340
00341
00342 pExtOrder = &Matrix->IntToExtRowMap[1];
00343 for (I = 1; I <= lSize; I++)
00344 { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
00345 { pElement = Matrix->FirstInRow[I];
00346
00347 while (pElement != NULL)
00348 { pElement->Real *= ScaleFactor;
00349 pElement = pElement->NextInRow;
00350 }
00351 }
00352 }
00353
00354
00355 pExtOrder = &Matrix->IntToExtColMap[1];
00356 for (I = 1; I <= lSize; I++)
00357 { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
00358 { pElement = Matrix->FirstInCol[I];
00359
00360 while (pElement != NULL)
00361 { pElement->Real *= ScaleFactor;
00362 pElement = pElement->NextInCol;
00363 }
00364 }
00365 }
00366 return;
00367
00368 #endif
00369 }
00370 #endif
00371
00372 #if spCOMPLEX && SCALING
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 static void ScaleComplexMatrix(MatrixPtr Matrix, double *RHS_ScaleFactors, double *SolutionScaleFactors )
00430 {
00431 ElementPtr pElement;
00432 int I, lSize, *pExtOrder;
00433 double ScaleFactor;
00434
00435
00436 lSize = Matrix->Size;
00437
00438
00439 #if !ARRAY_OFFSET
00440 --RHS_ScaleFactors;
00441 --SolutionScaleFactors;
00442 #endif
00443
00444
00445 pExtOrder = &Matrix->IntToExtRowMap[1];
00446 for (I = 1; I <= lSize; I++)
00447 { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
00448 { pElement = Matrix->FirstInRow[I];
00449
00450 while (pElement != NULL)
00451 { pElement->Real *= ScaleFactor;
00452 pElement->Imag *= ScaleFactor;
00453 pElement = pElement->NextInRow;
00454 }
00455 }
00456 }
00457
00458
00459 pExtOrder = &Matrix->IntToExtColMap[1];
00460 for (I = 1; I <= lSize; I++)
00461 { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
00462 { pElement = Matrix->FirstInCol[I];
00463
00464 while (pElement != NULL)
00465 { pElement->Real *= ScaleFactor;
00466 pElement->Imag *= ScaleFactor;
00467 pElement = pElement->NextInCol;
00468 }
00469 }
00470 }
00471 return;
00472 }
00473 #endif
00474
00475
00476
00477
00478
00479
00480
00481
00482 #if MULTIPLICATION
00483
00484
00485
00486
00487
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 void spMultiply(MatrixPtr Matrix, double* RHS, double* Solution)
00514 {
00515 ElementPtr pElement;
00516 double * Vector, Sum;
00517 int I, *pExtOrder;
00518
00519
00520 ASSERT( IS_SPARSE( Matrix ) && !Matrix->Factored );
00521 if (!Matrix->RowsLinked) spcLinkRows(Matrix);
00522
00523 #if spCOMPLEX
00524 if (Matrix->Complex)
00525 { ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
00526 return;
00527 }
00528 #endif
00529
00530 #if REAL
00531 #if !ARRAY_OFFSET
00532
00533 --RHS;
00534 --Solution;
00535 #endif
00536
00537
00538 Vector = Matrix->Intermediate;
00539 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
00540 for (I = Matrix->Size; I > 0; I--)
00541 Vector[I] = Solution[*(pExtOrder--)];
00542
00543 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
00544 for (I = Matrix->Size; I > 0; I--)
00545 {
00546 pElement = Matrix->FirstInRow[I];
00547 Sum = 0.0;
00548 while (pElement != NULL)
00549 {
00550 Sum += pElement->Real * Vector[pElement->Col];
00551 pElement = pElement->NextInRow;
00552 }
00553 RHS[*pExtOrder--] = Sum;
00554 }
00555 return;
00556 #endif
00557 }
00558 #endif
00559
00560
00561 #if MULTIPLICATION && TRANSPOSE
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592 void spMultTransposed(MatrixPtr Matrix, double *RHS, double *Solution)
00593 {
00594 ElementPtr pElement;
00595 double * Vector, Sum;
00596 int I, *pExtOrder;
00597
00598
00599
00600 ASSERT( IS_SPARSE( Matrix ) && !Matrix->Factored );
00601
00602 #if REAL
00603 #if !ARRAY_OFFSET
00604
00605 --RHS;
00606 --Solution;
00607 #endif
00608
00609
00610 Vector = Matrix->Intermediate;
00611 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
00612 for (I = Matrix->Size; I > 0; I--)
00613 Vector[I] = Solution[*(pExtOrder--)];
00614
00615 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
00616 for (I = Matrix->Size; I > 0; I--)
00617 { pElement = Matrix->FirstInCol[I];
00618 Sum = 0.0;
00619
00620 while (pElement != NULL)
00621 { Sum += pElement->Real * Vector[pElement->Row];
00622 pElement = pElement->NextInCol;
00623 }
00624 RHS[*pExtOrder--] = Sum;
00625 }
00626 return;
00627 #endif
00628 }
00629 #endif
00630
00631
00632
00633
00634
00635
00636
00637 #if spCOMPLEX && MULTIPLICATION && TRANSPOSE
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, double *RHS, double *Solution IMAG_VECTORS )
00673 {
00674 ElementPtr pElement;
00675 ComplexVector Vector;
00676 ComplexNumber Sum;
00677 int I, *pExtOrder;
00678
00679
00680
00681
00682 #if !ARRAY_OFFSET
00683 #if spSEPARATED_COMPLEX_VECTORS
00684 --RHS; --iRHS;
00685 --Solution; --iSolution;
00686 #else
00687 RHS -= 2; Solution -= 2;
00688 #endif
00689 #endif
00690
00691
00692 Vector = (ComplexVector)Matrix->Intermediate;
00693 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
00694
00695 #if spSEPARATED_COMPLEX_VECTORS
00696 for (I = Matrix->Size; I > 0; I--)
00697 { Vector[I].Real = Solution[*pExtOrder];
00698 Vector[I].Imag = iSolution[*(pExtOrder--)];
00699 }
00700 #else
00701 for (I = Matrix->Size; I > 0; I--)
00702 Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
00703 #endif
00704
00705 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
00706 for (I = Matrix->Size; I > 0; I--)
00707 { pElement = Matrix->FirstInCol[I];
00708 Sum.Real = Sum.Imag = 0.0;
00709
00710 while (pElement != NULL)
00711 {
00712 CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Row] );
00713 pElement = pElement->NextInCol;
00714 }
00715
00716 #if spSEPARATED_COMPLEX_VECTORS
00717 RHS[*pExtOrder] = Sum.Real;
00718 iRHS[*pExtOrder--] = Sum.Imag;
00719 #else
00720 ((ComplexVector)RHS)[*pExtOrder--] = Sum;
00721 #endif
00722 }
00723 return;
00724 }
00725 #endif
00726
00727 #if DETERMINANT
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768 #if spCOMPLEX
00769 void spDeterminant(MatrixPtr Matrix, int *pExponent, double *pDeterminant, double *piDeterminant )
00770 #else
00771 void spDeterminant(MatrixPtr Matrix, int *pExponent, double *pDeterminant )
00772 #endif
00773 {
00774 int I, Size;
00775
00776
00777
00778 #define NORM(a) (nr = ABSS((a).Real), ni = ABSS((a).Imag), MAX (nr,ni))
00779
00780
00781 ASSERT( IS_SPARSE( Matrix ) && IS_FACTORED(Matrix) );
00782 *pExponent = 0;
00783
00784 if (Matrix->Error == spSINGULAR)
00785 { *pDeterminant = 0.0;
00786 #if spCOMPLEX
00787 if (Matrix->Complex) *piDeterminant = 0.0;
00788 #endif
00789 return;
00790 }
00791
00792 Size = Matrix->Size;
00793 I = 0;
00794
00795 #if spCOMPLEX
00796 if (Matrix->Complex)
00797 { cDeterminant.Real = 1.0;
00798 cDeterminant.Imag = 0.0;
00799
00800 while (++I <= Size)
00801 { CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] );
00802 CMPLX_MULT_ASSIGN( cDeterminant, Pivot );
00803
00804
00805 Norm = NORM( cDeterminant );
00806 if (Norm != 0.0)
00807 { while (Norm >= 1.0e12)
00808 { cDeterminant.Real *= 1.0e-12;
00809 cDeterminant.Imag *= 1.0e-12;
00810 *pExponent += 12;
00811 Norm = NORM( cDeterminant );
00812 }
00813 while (Norm < 1.0e-12)
00814 { cDeterminant.Real *= 1.0e12;
00815 cDeterminant.Imag *= 1.0e12;
00816 *pExponent -= 12;
00817 Norm = NORM( cDeterminant );
00818 }
00819 }
00820 }
00821
00822
00823 Norm = NORM( cDeterminant );
00824 if (Norm != 0.0)
00825 { while (Norm >= 10.0)
00826 { cDeterminant.Real *= 0.1;
00827 cDeterminant.Imag *= 0.1;
00828 (*pExponent)++;
00829 Norm = NORM( cDeterminant );
00830 }
00831 while (Norm < 1.0)
00832 { cDeterminant.Real *= 10.0;
00833 cDeterminant.Imag *= 10.0;
00834 (*pExponent)--;
00835 Norm = NORM( cDeterminant );
00836 }
00837 }
00838 if (Matrix->NumberOfInterchangesIsOdd)
00839 CMPLX_NEGATE( cDeterminant );
00840
00841 *pDeterminant = cDeterminant.Real;
00842 *piDeterminant = cDeterminant.Imag;
00843 }
00844 #endif
00845 #if REAL && spCOMPLEX
00846 else
00847 #endif
00848 #if REAL
00849 {
00850 *pDeterminant = 1.0;
00851
00852 while (++I <= Size)
00853 { *pDeterminant /= Matrix->Diag[I]->Real;
00854
00855
00856 if (*pDeterminant != 0.0)
00857 { while (ABSS(*pDeterminant) >= 1.0e12)
00858 { *pDeterminant *= 1.0e-12;
00859 *pExponent += 12;
00860 }
00861 while (ABSS(*pDeterminant) < 1.0e-12)
00862 { *pDeterminant *= 1.0e12;
00863 *pExponent -= 12;
00864 }
00865 }
00866 }
00867
00868
00869 if (*pDeterminant != 0.0)
00870 { while (ABSS(*pDeterminant) >= 10.0)
00871 { *pDeterminant *= 0.1;
00872 (*pExponent)++;
00873 }
00874 while (ABSS(*pDeterminant) < 1.0)
00875 { *pDeterminant *= 10.0;
00876 (*pExponent)--;
00877 }
00878 }
00879 if (Matrix->NumberOfInterchangesIsOdd)
00880 *pDeterminant = -*pDeterminant;
00881 }
00882 #endif
00883 }
00884 #endif
00885
00886
00887
00888
00889
00890
00891
00892
00893 #if STRIP
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917 void spStripFills(MatrixPtr Matrix )
00918 {
00919 struct FillinListNodeStruct *pListNode;
00920
00921
00922 ASSERT( IS_SPARSE( Matrix ) );
00923 if (Matrix->Fillins == 0) return;
00924 Matrix->NeedsOrdering = 1;
00925 Matrix->Elements -= Matrix->Fillins;
00926 Matrix->Fillins = 0;
00927
00928
00929 { ElementPtr pFillin, pLastFillin;
00930
00931 pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
00932 Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
00933 Matrix->NextAvailFillin = pListNode->pFillinList;
00934
00935 while (pListNode != NULL)
00936 { pFillin = pListNode->pFillinList;
00937 pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]);
00938 while (pFillin <= pLastFillin)
00939 (pFillin++)->Row = 0;
00940 pListNode = pListNode->Next;
00941 }
00942 }
00943
00944
00945 { ElementPtr pElement, *ppElement;
00946 int I, Size = Matrix->Size;
00947
00948
00949 for (I = 1; I <= Size; I++)
00950 { ppElement = &(Matrix->FirstInCol[I]);
00951 while ((pElement = *ppElement) != NULL)
00952 { if (pElement->Row == 0)
00953 { *ppElement = pElement->NextInCol;
00954 if (Matrix->Diag[pElement->Col] == pElement)
00955 Matrix->Diag[pElement->Col] = NULL;
00956 }
00957 else
00958 ppElement = &pElement->NextInCol;
00959 }
00960 }
00961
00962
00963 for (I = 1; I <= Size; I++)
00964 { ppElement = &(Matrix->FirstInRow[I]);
00965 while ((pElement = *ppElement) != NULL)
00966 { if (pElement->Row == 0)
00967 *ppElement = pElement->NextInRow;
00968 else
00969 ppElement = &pElement->NextInRow;
00970 }
00971 }
00972 }
00973 return;
00974 }
00975 #endif
00976
00977 #if TRANSLATE && DELETE
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011 extern void spcColExchange(MatrixPtr Matrix, int Col1, int Col2 );
01012 extern void spcRowExchange(MatrixPtr Matrix, int Row1, int Row2 );
01013
01014 void spDeleteRowAndCol(MatrixPtr Matrix, int Row, int Col )
01015 {
01016 ElementPtr pElement, *ppElement, pLastElement, spcFindElementInCol();
01017 int Size, ExtRow, ExtCol;
01018
01019 ASSERT( IS_SPARSE(Matrix) && Row > 0 && Col > 0 );
01020 ASSERT( Row <= Matrix->ExtSize && Col <= Matrix->ExtSize );
01021
01022 Size = Matrix->Size;
01023 ExtRow = Row;
01024 ExtCol = Col;
01025 if (!Matrix->RowsLinked) spcLinkRows( Matrix );
01026
01027 Row = Matrix->ExtToIntRowMap[Row];
01028 Col = Matrix->ExtToIntColMap[Col];
01029 ASSERT( Row > 0 && Col > 0 );
01030
01031
01032 if (Row != Size) spcRowExchange( Matrix, Row, Size );
01033
01034
01035 if (Col != Size) spcColExchange( Matrix, Col, Size );
01036
01037
01038 if (Row == Col)
01039 SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] )
01040 else
01041 {
01042 Matrix->Diag[Row] = spcFindElementInCol(Matrix, Matrix->FirstInCol+Row, Row, Row, 0 );
01043 Matrix->Diag[Col] = spcFindElementInCol(Matrix, Matrix->FirstInCol+Col, Col, Col, 0);
01044 }
01045
01046
01047
01048
01049
01050 pLastElement = Matrix->FirstInRow[ Size ];
01051 while (pLastElement != NULL)
01052 { ppElement = &(Matrix->FirstInCol[ pLastElement->Col ]);
01053 while ((pElement = *ppElement) != NULL)
01054 { if (pElement == pLastElement)
01055 *ppElement = NULL;
01056 else
01057 ppElement = &pElement->NextInCol;
01058 }
01059 pLastElement = pLastElement->NextInRow;
01060 }
01061
01062
01063 pLastElement = Matrix->FirstInCol[ Size ];
01064 while (pLastElement != NULL)
01065 { ppElement = &(Matrix->FirstInRow[ pLastElement->Row ]);
01066 while ((pElement = *ppElement) != NULL)
01067 { if (pElement == pLastElement)
01068 *ppElement = NULL;
01069 else
01070 ppElement = &pElement->NextInRow;
01071 }
01072 pLastElement = pLastElement->NextInCol;
01073 }
01074
01075
01076 Matrix->Size = Size - 1;
01077 Matrix->Diag[Size] = NULL;
01078 Matrix->FirstInRow[Size] = NULL;
01079 Matrix->FirstInCol[Size] = NULL;
01080 Matrix->CurrentSize--;
01081 Matrix->ExtToIntRowMap[ExtRow] = -1;
01082 Matrix->ExtToIntColMap[ExtCol] = -1;
01083 Matrix->NeedsOrdering = 1;
01084
01085 return;
01086 }
01087 #endif
01088
01089 #if PSEUDOCONDITION
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114 double spPseudoCondition(MatrixPtr Matrix)
01115 {
01116 int I;
01117 ArrayOfElementPtrs Diag;
01118 double MaxPivot, MinPivot, Mag;
01119
01120
01121
01122 ASSERT( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
01123 if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG)
01124 return 0.0;
01125
01126 Diag = Matrix->Diag;
01127 MaxPivot = MinPivot = ELEMENT_MAG( Diag[1] );
01128 for (I = 2; I <= Matrix->Size; I++)
01129 { Mag = ELEMENT_MAG( Diag[I] );
01130 if (Mag > MaxPivot)
01131 MaxPivot = Mag;
01132 else if (Mag < MinPivot)
01133 MinPivot = Mag;
01134 }
01135 ASSERT( MaxPivot > 0.0);
01136 return MaxPivot / MinPivot;
01137 }
01138 #endif
01139
01140 #if CONDITION
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193 double spCondition(MatrixPtr Matrix, double NormOfMatrix, int *pError)
01194 {
01195 ElementPtr pElement, pPivot;
01196 double * T, *Tm;
01197 int I, K, Row;
01198 int Size;
01199 double E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor,
01200 Linpack, OLeary, InvNormOfInverse;
01201 #define SLACK 1e4
01202
01203
01204
01205 ASSERT( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
01206 *pError = Matrix->Error;
01207 if (Matrix->Error >= spFATAL) return 0.0;
01208 if (NormOfMatrix == 0.0)
01209 { *pError = spSINGULAR;
01210 return 0.0;
01211 }
01212
01213 #if spCOMPLEX
01214 if (Matrix->Complex)
01215 return ComplexCondition( Matrix, NormOfMatrix, pError );
01216 #endif
01217
01218 #if REAL
01219 Size = Matrix->Size;
01220 T = Matrix->Intermediate;
01221 #if spCOMPLEX
01222 Tm = Matrix->Intermediate + Size;
01223 #else
01224 Tm = ALLOC( double, Size+1 );
01225 if (Tm == NULL)
01226 {
01227 *pError = spNO_MEMORY;
01228 return 0.0;
01229 }
01230 #endif
01231 for (I = Size; I > 0; I--) T[I] = 0.0;
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241 E = 1.0;
01242 for (I = 1; I <= Size; I++)
01243 {
01244 pPivot = Matrix->Diag[I];
01245 if (T[I] < 0.0) Em = -E; else Em = E;
01246 Wm = (Em + T[I]) * pPivot->Real;
01247 if (ABSS(Wm) > SLACK)
01248 {
01249 ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABSS(Wm) );
01250 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
01251 E *= ScaleFactor;
01252 Em *= ScaleFactor;
01253 Wm = (Em + T[I]) * pPivot->Real;
01254 }
01255 Wp = (T[I] - Em) * pPivot->Real;
01256 ASp = ABSS(T[I] - Em);
01257 ASm = ABSS(Em + T[I]);
01258
01259
01260 pElement = pPivot->NextInCol;
01261 while (pElement != NULL)
01262 {
01263 Row = pElement->Row;
01264 Tm[Row] = T[Row] - (Wm * pElement->Real);
01265 T[Row] -= (Wp * pElement->Real);
01266 ASp += ABSS(T[Row]);
01267 ASm += ABSS(Tm[Row]);
01268 pElement = pElement->NextInCol;
01269 }
01270
01271
01272 if (ASm > ASp)
01273 {
01274 T[I] = Wm;
01275 pElement = pPivot->NextInCol;
01276 while (pElement != NULL)
01277 {
01278 T[pElement->Row] = Tm[pElement->Row];
01279 pElement = pElement->NextInCol;
01280 }
01281 }
01282 else T[I] = Wp;
01283 }
01284
01285
01286 for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABSS(T[I]);
01287 ScaleFactor = 1.0 / (SLACK * ASw);
01288 if (ScaleFactor < 0.5)
01289 {
01290 for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
01291 E *= ScaleFactor;
01292 }
01293
01294
01295 for (I = Size; I >= 1; I--)
01296 { pElement = Matrix->Diag[I]->NextInRow;
01297 while (pElement != NULL)
01298 {
01299 T[I] -= pElement->Real * T[pElement->Col];
01300 pElement = pElement->NextInRow;
01301 }
01302 if (ABSS(T[I]) > SLACK)
01303 {
01304 ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABSS(T[I]) );
01305 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
01306 E *= ScaleFactor;
01307 }
01308 }
01309
01310
01311 for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABSS(T[I]);
01312 ScaleFactor = 1.0 / (SLACK * ASy);
01313 if (ScaleFactor < 0.5)
01314 {
01315 for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
01316 ASy = 1.0 / SLACK;
01317 E *= ScaleFactor;
01318 }
01319
01320
01321 for (MaxY = 0.0, I = Size; I > 0; I--)
01322 if (MaxY < ABSS(T[I])) MaxY = ABSS(T[I]);
01323
01324
01325
01326
01327
01328
01329
01330 for (I = 1; I <= Size; I++)
01331 {
01332 pElement = Matrix->Diag[I]->NextInRow;
01333 while (pElement != NULL)
01334 {
01335 T[pElement->Col] -= T[I] * pElement->Real;
01336 pElement = pElement->NextInRow;
01337 }
01338 if (ABSS(T[I]) > SLACK)
01339 {
01340 ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABSS(T[I]) );
01341 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
01342 ASy *= ScaleFactor;
01343 }
01344 }
01345
01346
01347 for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABSS(T[I]);
01348 ScaleFactor = 1.0 / (SLACK * ASv);
01349 if (ScaleFactor < 0.5)
01350 {
01351 for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
01352 ASy *= ScaleFactor;
01353 }
01354
01355
01356 for (I = Size; I >= 1; I--)
01357 {
01358 pPivot = Matrix->Diag[I];
01359 pElement = pPivot->NextInCol;
01360 while (pElement != NULL)
01361 {
01362 T[I] -= pElement->Real * T[pElement->Row];
01363 pElement = pElement->NextInCol;
01364 }
01365 T[I] *= pPivot->Real;
01366 if (ABSS(T[I]) > SLACK)
01367 {
01368 ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABSS(T[I]) );
01369 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
01370 ASy *= ScaleFactor;
01371 }
01372 }
01373
01374
01375 for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABSS(T[I]);
01376
01377 #if !spCOMPLEX
01378 FREE( Tm );
01379 #endif
01380
01381 Linpack = ASy / ASz;
01382 OLeary = E / MaxY;
01383 InvNormOfInverse = MIN( Linpack, OLeary );
01384 return InvNormOfInverse / NormOfMatrix;
01385 #endif
01386 }
01387
01388 #if spCOMPLEX
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410 static double ComplexCondition(MatrixPtr Matrix, double NormOfMatrix, int *pError)
01411 {
01412 ElementPtr pElement;
01413 ComplexVector T, Tm;
01414 int I, K, Row, Size;
01415 ElementPtr pPivot;
01416 double E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor, Linpack, OLeary, InvNormOfInverse;
01417 ComplexNumber Wp, Wm;
01418
01419 Size = Matrix->Size;
01420 T = (ComplexVector)Matrix->Intermediate;
01421 Tm = ALLOC( ComplexNumber, Size+1 );
01422 if (Tm == NULL)
01423 {
01424 *pError = spNO_MEMORY;
01425 return 0.0;
01426 }
01427 for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0;
01428
01429
01430
01431
01432
01433
01434
01435 E = 1.0;
01436 for (I = 1; I <= Size; I++)
01437 {
01438 pPivot = Matrix->Diag[I];
01439 if (T[I].Real < 0.0) Em = -E; else Em = E;
01440 Wm = T[I];
01441 Wm.Real += Em;
01442 ASm = CMPLX_1_NORM( Wm );
01443 CMPLX_MULT_ASSIGN( Wm, *pPivot );
01444 if (CMPLX_1_NORM(Wm) > SLACK)
01445 {
01446 ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) );
01447 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
01448 E *= ScaleFactor;
01449 Em *= ScaleFactor;
01450 ASm *= ScaleFactor;
01451 SCLR_MULT_ASSIGN( Wm, ScaleFactor );
01452 }
01453 Wp = T[I];
01454 Wp.Real -= Em;
01455 ASp = CMPLX_1_NORM( Wp );
01456 CMPLX_MULT_ASSIGN( Wp, *pPivot );
01457
01458
01459 pElement = pPivot->NextInCol;
01460 while (pElement != NULL)
01461 {
01462 Row = pElement->Row;
01463
01464 CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] );
01465
01466 CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement );
01467 ASp += CMPLX_1_NORM(T[Row]);
01468 ASm += CMPLX_1_NORM(Tm[Row]);
01469 pElement = pElement->NextInCol;
01470 }
01471
01472
01473 if (ASm > ASp)
01474 {
01475 T[I] = Wm;
01476 pElement = pPivot->NextInCol;
01477 while (pElement != NULL)
01478 {
01479 T[pElement->Row] = Tm[pElement->Row];
01480 pElement = pElement->NextInCol;
01481 }
01482 }
01483 else T[I] = Wp;
01484 }
01485
01486
01487 for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]);
01488 ScaleFactor = 1.0 / (SLACK * ASw);
01489 if (ScaleFactor < 0.5)
01490 {
01491 for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
01492 E *= ScaleFactor;
01493 }
01494
01495
01496 for (I = Size; I >= 1; I--)
01497 {
01498 pElement = Matrix->Diag[I]->NextInRow;
01499 while (pElement != NULL)
01500 {
01501 CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement );
01502 pElement = pElement->NextInRow;
01503 }
01504 if (CMPLX_1_NORM(T[I]) > SLACK)
01505 {
01506 ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
01507 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
01508 E *= ScaleFactor;
01509 }
01510 }
01511
01512
01513 for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]);
01514 ScaleFactor = 1.0 / (SLACK * ASy);
01515 if (ScaleFactor < 0.5)
01516 {
01517 for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
01518 ASy = 1.0 / SLACK;
01519 E *= ScaleFactor;
01520 }
01521
01522
01523 for (MaxY = 0.0, I = Size; I > 0; I--)
01524 if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]);
01525
01526
01527
01528
01529
01530
01531
01532 for (I = 1; I <= Size; I++)
01533 {
01534 pElement = Matrix->Diag[I]->NextInRow;
01535 while (pElement != NULL)
01536 {
01537 CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement );
01538 pElement = pElement->NextInRow;
01539 }
01540 if (CMPLX_1_NORM(T[I]) > SLACK)
01541 {
01542 ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
01543 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
01544 ASy *= ScaleFactor;
01545 }
01546 }
01547
01548
01549 for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]);
01550 ScaleFactor = 1.0 / (SLACK * ASv);
01551 if (ScaleFactor < 0.5)
01552 {
01553 for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
01554 ASy *= ScaleFactor;
01555 }
01556
01557
01558 for (I = Size; I >= 1; I--)
01559 {
01560 pPivot = Matrix->Diag[I];
01561 pElement = pPivot->NextInCol;
01562 while (pElement != NULL)
01563 {
01564 CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement );
01565 pElement = pElement->NextInCol;
01566 }
01567 CMPLX_MULT_ASSIGN( T[I], *pPivot );
01568 if (CMPLX_1_NORM(T[I]) > SLACK)
01569 {
01570 ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
01571 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
01572 ASy *= ScaleFactor;
01573 }
01574 }
01575
01576
01577 for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]);
01578
01579 FREE( Tm );
01580
01581 Linpack = ASy / ASz;
01582 OLeary = E / MaxY;
01583 InvNormOfInverse = MIN( Linpack, OLeary );
01584 return InvNormOfInverse / NormOfMatrix;
01585 }
01586 #endif
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604 double spNorm(MatrixPtr Matrix)
01605 {
01606 ElementPtr pElement;
01607 int I;
01608 double Max = 0.0, AbsRowSum;
01609
01610 ASSERT( IS_SPARSE(Matrix) && !IS_FACTORED(Matrix) );
01611 if (!Matrix->RowsLinked) spcLinkRows( Matrix );
01612
01613
01614 #if REAL
01615 if (!Matrix->Complex)
01616 {
01617 for (I = Matrix->Size; I > 0; I--)
01618 {
01619 pElement = Matrix->FirstInRow[I];
01620 AbsRowSum = 0.0;
01621 while (pElement != NULL)
01622 {
01623 AbsRowSum += ABSS( pElement->Real );
01624 pElement = pElement->NextInRow;
01625 }
01626 if (Max < AbsRowSum) Max = AbsRowSum;
01627 }
01628 }
01629 #endif
01630 return Max;
01631 }
01632 #endif
01633
01634 #if STABILITY
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701 double spLargestElement(MatrixPtr Matrix)
01702 {
01703 int I;
01704 double Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0;
01705 double Pivot;
01706
01707 ElementPtr pElement, pDiag;
01708
01709
01710 ASSERT( IS_SPARSE(Matrix) );
01711
01712 #if REAL
01713 if (Matrix->Factored && !Matrix->Complex)
01714 { if (Matrix->Error == spSINGULAR) return 0.0;
01715
01716
01717 for (I = 1; I <= Matrix->Size; I++)
01718 { pDiag = Matrix->Diag[I];
01719
01720
01721 Pivot = 1.0 / pDiag->Real;
01722 Mag = ABSS( Pivot );
01723 if (Mag > MaxRow) MaxRow = Mag;
01724 pElement = Matrix->FirstInRow[I];
01725 while (pElement != pDiag)
01726 { Mag = ABSS( pElement->Real );
01727 if (Mag > MaxRow) MaxRow = Mag;
01728 pElement = pElement->NextInRow;
01729 }
01730
01731
01732 pElement = Matrix->FirstInCol[I];
01733 AbsColSum = 1.0;
01734 while (pElement != pDiag)
01735 { AbsColSum += ABSS( pElement->Real );
01736 pElement = pElement->NextInCol;
01737 }
01738 if (AbsColSum > MaxCol) MaxCol = AbsColSum;
01739 }
01740 }
01741 else if (!Matrix->Complex)
01742 { for (I = 1; I <= Matrix->Size; I++)
01743 { pElement = Matrix->FirstInCol[I];
01744 while (pElement != NULL)
01745 { Mag = ABSS( pElement->Real );
01746 if (Mag > Max) Max = Mag;
01747 pElement = pElement->NextInCol;
01748 }
01749 }
01750 return Max;
01751 }
01752 #endif
01753 #if spCOMPLEX
01754 if (Matrix->Factored && Matrix->Complex)
01755 { if (Matrix->Error == spSINGULAR) return 0.0;
01756
01757
01758 for (I = 1; I <= Matrix->Size; I++)
01759 { pDiag = Matrix->Diag[I];
01760
01761
01762 CMPLX_RECIPROCAL( cPivot, *pDiag );
01763 Mag = CMPLX_1_NORM( cPivot );
01764 if (Mag > MaxRow) MaxRow = Mag;
01765 pElement = Matrix->FirstInRow[I];
01766 while (pElement != pDiag)
01767 { Mag = CMPLX_1_NORM( *pElement );
01768 if (Mag > MaxRow) MaxRow = Mag;
01769 pElement = pElement->NextInRow;
01770 }
01771
01772
01773 pElement = Matrix->FirstInCol[I];
01774 AbsColSum = 1.0;
01775 while (pElement != pDiag)
01776 { AbsColSum += CMPLX_1_NORM( *pElement );
01777 pElement = pElement->NextInCol;
01778 }
01779 if (AbsColSum > MaxCol) MaxCol = AbsColSum;
01780 }
01781 }
01782 else if (Matrix->Complex)
01783 { for (I = 1; I <= Matrix->Size; I++)
01784 { pElement = Matrix->FirstInCol[I];
01785 while (pElement != NULL)
01786 { Mag = CMPLX_1_NORM( *pElement );
01787 if (Mag > Max) Max = Mag;
01788 pElement = pElement->NextInCol;
01789 }
01790 }
01791 return Max;
01792 }
01793 #endif
01794 return MaxRow * MaxCol;
01795 }
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812 double spRoundoff(MatrixPtr Matrix, double Rho)
01813 {
01814 ElementPtr pElement;
01815 int Count, I, MaxCount = 0;
01816 double Reid, Gear;
01817
01818 ASSERT( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
01819
01820
01821 if (Rho < 0.0) Rho = spLargestElement(Matrix);
01822
01823
01824 if (Matrix->MaxRowCountInLowerTri < 0)
01825 { for (I = Matrix->Size; I > 0; I--)
01826 { pElement = Matrix->FirstInRow[I];
01827 Count = 0;
01828 while (pElement->Col < I)
01829 { Count++;
01830 pElement = pElement->NextInRow;
01831 }
01832 if (Count > MaxCount) MaxCount = Count;
01833 }
01834 Matrix->MaxRowCountInLowerTri = MaxCount;
01835 }
01836 else MaxCount = Matrix->MaxRowCountInLowerTri;
01837
01838
01839 Gear = 1.01*((MaxCount + 1) * Matrix->RelThreshold + 1.0) * SQR(MaxCount);
01840 Reid = 3.01 * Matrix->Size;
01841
01842 if (Gear < Reid)
01843 return (DBL_EPSILON * Rho * Gear);
01844 else
01845 return (DBL_EPSILON * Rho * Reid);
01846 }
01847 #endif