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 #define spINSIDE_SPARSE
00060 #include "spconfig.h"
00061 #include "spmatrix.h"
00062 extern ElementPtr spcCreateElement(MatrixPtr Matrix, int Row, int Col, ElementPtr *LastAddr, BOOLEAN Fillin);
00063 static ElementPtr SearchDiagonal(MatrixPtr Matrix, int Step);
00064 static ElementPtr QuicklySearchDiagonal(MatrixPtr Matrix, int Step);
00065 static ElementPtr SearchEntireMatrix(MatrixPtr Matrix, int Step);
00066 static double FindLargestInCol(ElementPtr pElement);
00067 static ElementPtr CreateFillin(MatrixPtr Matrix, int Row, int Col);
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
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 void spcColExchange(MatrixPtr Matrix, int Col1, int Col2 );
00162 static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot);
00163 static void ExchangeRowsAndCols(MatrixPtr Matrix, ElementPtr pPivot, int Step );
00164 static void CountMarkowitz(MatrixPtr Matrix, double *RHS, int Step );
00165 static void CreateInternalVectors(MatrixPtr Matrix);
00166 static void MarkowitzProducts(MatrixPtr Matrix, int Step);
00167 static int MatrixIsSingular(MatrixPtr Matrix, int Step);
00168 static void UpdateMarkowitzNumbers(MatrixPtr Matrix, ElementPtr pPivot);
00169 extern void spcLinkRows(MatrixPtr Matrix);
00170 static int ZeroPivot(MatrixPtr Matrix, int Step );
00171 static double FindBiggestInColExclude(MatrixPtr Matrix, ElementPtr pElement, int Step);
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 static ElementPtr SearchForSingleton(MatrixPtr Matrix, int Step)
00209 {
00210 ElementPtr ChosenPivot;
00211 int I;
00212 long *pMarkowitzProduct;
00213 int Singletons;
00214 double PivotMag;
00215
00216
00217 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+1]);
00218 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
00219
00220
00221
00222 Singletons = Matrix->Singletons--;
00223
00224
00225
00226
00227
00228
00229 Matrix->MarkowitzProd[Step-1] = 0;
00230
00231 while (Singletons-- > 0)
00232 {
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 while ( *pMarkowitzProduct-- )
00255 {
00256
00257
00258
00259
00260
00261 }
00262 I = (int)(pMarkowitzProduct - Matrix->MarkowitzProd + 1);
00263
00264
00265 if (I < Step) break;
00266 if (I > Matrix->Size) I = Step;
00267
00268
00269 if ((ChosenPivot = Matrix->Diag[I]) != NULL)
00270 {
00271
00272 PivotMag = ELEMENT_MAG(ChosenPivot);
00273 if
00274 ( PivotMag > Matrix->AbsThreshold &&
00275 PivotMag > Matrix->RelThreshold *
00276 FindBiggestInColExclude( Matrix, ChosenPivot, Step )
00277 ) return ChosenPivot;
00278 }
00279 else
00280 {
00281
00282 if (Matrix->MarkowitzCol[I] == 0)
00283 { ChosenPivot = Matrix->FirstInCol[I];
00284 while ((ChosenPivot != NULL) && (ChosenPivot->Row < Step))
00285 ChosenPivot = ChosenPivot->NextInCol;
00286 PivotMag = ELEMENT_MAG(ChosenPivot);
00287 if
00288 ( PivotMag > Matrix->AbsThreshold &&
00289 PivotMag > Matrix->RelThreshold *
00290 FindBiggestInColExclude( Matrix, ChosenPivot,
00291 Step )
00292 ) return ChosenPivot;
00293 else
00294 { if (Matrix->MarkowitzRow[I] == 0)
00295 { ChosenPivot = Matrix->FirstInRow[I];
00296 while((ChosenPivot != NULL) && (ChosenPivot->Col<Step))
00297 ChosenPivot = ChosenPivot->NextInRow;
00298 PivotMag = ELEMENT_MAG(ChosenPivot);
00299 if
00300 ( PivotMag > Matrix->AbsThreshold &&
00301 PivotMag > Matrix->RelThreshold *
00302 FindBiggestInColExclude( Matrix,
00303 ChosenPivot,
00304 Step )
00305 ) return ChosenPivot;
00306 }
00307 }
00308 }
00309 else
00310 { ChosenPivot = Matrix->FirstInRow[I];
00311 while ((ChosenPivot != NULL) && (ChosenPivot->Col < Step))
00312 ChosenPivot = ChosenPivot->NextInRow;
00313 PivotMag = ELEMENT_MAG(ChosenPivot);
00314 if
00315 ( PivotMag > Matrix->AbsThreshold &&
00316 PivotMag > Matrix->RelThreshold *
00317 FindBiggestInColExclude( Matrix, ChosenPivot,
00318 Step )
00319 ) return ChosenPivot;
00320 }
00321 }
00322
00323 }
00324
00325
00326
00327
00328
00329 Matrix->Singletons++;
00330 return NULL;
00331 }
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 static ElementPtr SearchForPivot(MatrixPtr Matrix, int Step, int DiagPivoting)
00373 {
00374 ElementPtr ChosenPivot;
00375
00376
00377 if (Matrix->Singletons)
00378 { ChosenPivot = SearchForSingleton( Matrix, Step );
00379 if (ChosenPivot != NULL)
00380 { Matrix->PivotSelectionMethod = 's';
00381 return ChosenPivot;
00382 }
00383 }
00384
00385 #if DIAGONAL_PIVOTING
00386 if (DiagPivoting)
00387 {
00388
00389
00390
00391
00392
00393
00394 ChosenPivot = QuicklySearchDiagonal( Matrix, Step );
00395 if (ChosenPivot != NULL)
00396 {
00397 Matrix->PivotSelectionMethod = 'q';
00398 return ChosenPivot;
00399 }
00400
00401
00402
00403
00404
00405 ChosenPivot = SearchDiagonal( Matrix, Step );
00406 if (ChosenPivot != NULL)
00407 { Matrix->PivotSelectionMethod = 'd';
00408 return ChosenPivot;
00409 }
00410 }
00411 #endif
00412
00413
00414 ChosenPivot = SearchEntireMatrix( Matrix, Step );
00415 Matrix->PivotSelectionMethod = 'e';
00416
00417 return ChosenPivot;
00418 }
00419
00420 int spOrderAndFactor(MatrixPtr Matrix, double *RHS, double RelThreshold, double AbsThreshold, BOOLEAN DiagPivoting)
00421 {
00422 ElementPtr pPivot;
00423 int Step, Size, ReorderingRequired;
00424 double LargestInCol;
00425
00426
00427 ASSERT( IS_VALID(Matrix) && !Matrix->Factored);
00428
00429 Matrix->Error = spOKAY;
00430 Size = Matrix->Size;
00431 if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold;
00432 if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold;
00433 Matrix->RelThreshold = RelThreshold;
00434 if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold;
00435 Matrix->AbsThreshold = AbsThreshold;
00436 ReorderingRequired = 0;
00437
00438 if (!Matrix->NeedsOrdering)
00439 {
00440
00441 for (Step = 1; Step <= Size; Step++)
00442 { pPivot = Matrix->Diag[Step];
00443 LargestInCol = FindLargestInCol(pPivot->NextInCol);
00444 if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
00445 {
00446 RealRowColElimination( Matrix, pPivot );
00447 }
00448 else
00449 { ReorderingRequired = 1;
00450 break;
00451 }
00452 }
00453 if (!ReorderingRequired)
00454 goto Done;
00455 else
00456 {
00457
00458
00459
00460
00461
00462 #if ANNOTATE >= ON_STRANGE_BEHAVIOR
00463 printf("Reordering, Step = %1d\n", Step);
00464 #endif
00465 }
00466 }
00467 else
00468 {
00469
00470
00471
00472
00473
00474
00475 Step = 1;
00476 if (!Matrix->RowsLinked)
00477 spcLinkRows( Matrix );
00478 if (!Matrix->InternalVectorsAllocated)
00479 CreateInternalVectors( Matrix );
00480 if (Matrix->Error >= spFATAL)
00481 return Matrix->Error;
00482 }
00483
00484
00485 CountMarkowitz( Matrix, RHS, Step );
00486 MarkowitzProducts( Matrix, Step );
00487 Matrix->MaxRowCountInLowerTri = -1;
00488
00489
00490 for (; Step <= Size; Step++)
00491 { pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
00492 if (pPivot == NULL) return MatrixIsSingular( Matrix, Step );
00493 ExchangeRowsAndCols( Matrix, pPivot, Step );
00494
00495 RealRowColElimination( Matrix, pPivot );
00496
00497 if (Matrix->Error >= spFATAL) return Matrix->Error;
00498 UpdateMarkowitzNumbers( Matrix, pPivot );
00499
00500 #if ANNOTATE == FULL
00501 WriteStatus( Matrix, Step );
00502 #endif
00503 }
00504
00505 Done:
00506 Matrix->NeedsOrdering = 0;
00507 Matrix->Reordered = 1;
00508 Matrix->Factored = 1;
00509
00510 return Matrix->Error;
00511 }
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543 int spFactor(MatrixPtr Matrix)
00544 {
00545 ElementPtr pElement, pColumn;
00546 int Step, Size;
00547 double Mult;
00548
00549 ASSERT( IS_VALID(Matrix) && !Matrix->Factored);
00550
00551 if (Matrix->NeedsOrdering)
00552 return spOrderAndFactor(Matrix, (double*)NULL, 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT);
00553 if (!Matrix->Partitioned) spPartition(Matrix, spDEFAULT_PARTITION );
00554
00555 #if REAL
00556 Size = Matrix->Size;
00557
00558 if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
00559 Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;
00560
00561
00562 for (Step = 2; Step <= Size; Step++)
00563 { if (Matrix->DoRealDirect[Step])
00564 {
00565 double *Dest = (double *)Matrix->Intermediate;
00566
00567
00568 pElement = Matrix->FirstInCol[Step];
00569 while (pElement != NULL)
00570 { Dest[pElement->Row] = pElement->Real;
00571 pElement = pElement->NextInCol;
00572 }
00573
00574
00575 pColumn = Matrix->FirstInCol[Step];
00576 while (pColumn->Row < Step)
00577 { pElement = Matrix->Diag[pColumn->Row];
00578 pColumn->Real = Dest[pColumn->Row] * pElement->Real;
00579 while ((pElement = pElement->NextInCol) != NULL)
00580 Dest[pElement->Row] -= pColumn->Real * pElement->Real;
00581 pColumn = pColumn->NextInCol;
00582 }
00583
00584
00585 pElement = Matrix->Diag[Step]->NextInCol;
00586 while (pElement != NULL)
00587 { pElement->Real = Dest[pElement->Row];
00588 pElement = pElement->NextInCol;
00589 }
00590
00591
00592 if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step );
00593 Matrix->Diag[Step]->Real = 1.0 / Dest[Step];
00594 }
00595 else
00596 {
00597 double **pDest = (double **)Matrix->Intermediate;
00598
00599
00600 pElement = Matrix->FirstInCol[Step];
00601 while (pElement != NULL)
00602 { pDest[pElement->Row] = &pElement->Real;
00603 pElement = pElement->NextInCol;
00604 }
00605
00606
00607 pColumn = Matrix->FirstInCol[Step];
00608 while (pColumn->Row < Step)
00609 { pElement = Matrix->Diag[pColumn->Row];
00610 Mult = (*pDest[pColumn->Row] *= pElement->Real);
00611 while ((pElement = pElement->NextInCol) != NULL)
00612 *pDest[pElement->Row] -= Mult * pElement->Real;
00613 pColumn = pColumn->NextInCol;
00614 }
00615
00616
00617 if (Matrix->Diag[Step]->Real == 0.0)
00618 return ZeroPivot( Matrix, Step );
00619 Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real;
00620 }
00621 }
00622
00623 Matrix->Factored = 1;
00624 return (Matrix->Error = spOKAY);
00625 #endif
00626 }
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
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 void spPartition(MatrixPtr Matrix, int Mode)
00666 {
00667 ElementPtr pElement, pColumn;
00668 int Step, Size, *Nc, *No, *Nm;
00669 BOOLEAN *DoRealDirect, *DoCmplxDirect;
00670
00671 ASSERT( IS_SPARSE( Matrix ) );
00672 if (Matrix->Partitioned) return;
00673 Size = Matrix->Size;
00674 DoRealDirect = Matrix->DoRealDirect;
00675 DoCmplxDirect = Matrix->DoCmplxDirect;
00676 Matrix->Partitioned = 1;
00677
00678
00679 if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION;
00680 if (Mode == spDIRECT_PARTITION)
00681 {
00682 for (Step = 1; Step <= Size; Step++)
00683 DoRealDirect[Step] = 1;
00684 return;
00685 }
00686 else if (Mode == spINDIRECT_PARTITION)
00687 { for (Step = 1; Step <= Size; Step++)
00688 DoRealDirect[Step] = 0;
00689 return;
00690 }
00691 else ASSERT( Mode == spAUTO_PARTITION );
00692
00693
00694 Nc = (int *)Matrix->MarkowitzRow;
00695 No = (int *)Matrix->MarkowitzCol;
00696 Nm = (int *)Matrix->MarkowitzProd;
00697
00698
00699 for (Step = 1; Step <= Size; Step++)
00700 { Nc[Step] = No[Step] = Nm[Step] = 0;
00701
00702 pElement = Matrix->FirstInCol[Step];
00703 while (pElement != NULL)
00704 { Nc[Step]++;
00705 pElement = pElement->NextInCol;
00706 }
00707
00708 pColumn = Matrix->FirstInCol[Step];
00709 while (pColumn->Row < Step)
00710 { pElement = Matrix->Diag[pColumn->Row];
00711 Nm[Step]++;
00712 while ((pElement = pElement->NextInCol) != NULL)
00713 No[Step]++;
00714 pColumn = pColumn->NextInCol;
00715 }
00716 }
00717
00718 for (Step = 1; Step <= Size; Step++)
00719 {
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732 #define generic
00733 #ifdef hp9000s300
00734 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
00735 #undef generic
00736 #endif
00737
00738 #ifdef vax
00739 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
00740 #undef generic
00741 #endif
00742
00743 #ifdef generic
00744 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
00745 #undef generic
00746 #endif
00747 }
00748
00749 #if (ANNOTATE == FULL)
00750 { int Ops = 0;
00751 for (Step = 1; Step <= Size; Step++)
00752 Ops += No[Step];
00753 printf("Operation count for inner loop of factorization = %d.\n", Ops);
00754 }
00755 #endif
00756 return;
00757 }
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 static void CreateInternalVectors(MatrixPtr Matrix)
00777 {
00778 int Size;
00779
00780
00781
00782 Size= Matrix->Size;
00783
00784 if (Matrix->MarkowitzRow == NULL)
00785 { if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL)
00786 Matrix->Error = spNO_MEMORY;
00787 }
00788 if (Matrix->MarkowitzCol == NULL)
00789 { if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL)
00790 Matrix->Error = spNO_MEMORY;
00791 }
00792 if (Matrix->MarkowitzProd == NULL)
00793 { if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL)
00794 Matrix->Error = spNO_MEMORY;
00795 }
00796
00797
00798 #if REAL
00799 if (Matrix->DoRealDirect == NULL)
00800 { if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
00801 Matrix->Error = spNO_MEMORY;
00802 }
00803 #endif
00804 #if spCOMPLEX
00805 if (Matrix->DoCmplxDirect == NULL)
00806 { if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
00807 Matrix->Error = spNO_MEMORY;
00808 }
00809 #endif
00810
00811
00812 #if spCOMPLEX
00813 if (Matrix->Intermediate == NULL)
00814 { if ((Matrix->Intermediate = ALLOC(double,2*(Size+1))) == NULL)
00815 Matrix->Error = spNO_MEMORY;
00816 }
00817 #else
00818 if (Matrix->Intermediate == NULL)
00819 { if ((Matrix->Intermediate = ALLOC(double, Size+1)) == NULL)
00820 Matrix->Error = spNO_MEMORY;
00821 }
00822 #endif
00823
00824 if (Matrix->Error != spNO_MEMORY)
00825 Matrix->InternalVectorsAllocated = 1;
00826 return;
00827 }
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862 static void CountMarkowitz(MatrixPtr Matrix, double *RHS, int Step )
00863 {
00864 int Count, I, Size = Matrix->Size;
00865 ElementPtr pElement;
00866 int ExtRow;
00867
00868
00869
00870
00871 #if !ARRAY_OFFSET
00872 #if spSEPARATED_COMPLEX_VECTORS || !spCOMPLEX
00873 if (RHS != NULL) --RHS;
00874 #else
00875 if (RHS != NULL)
00876 { if (Matrix->Complex) RHS -= 2;
00877 else --RHS;
00878 }
00879 #endif
00880 #endif
00881
00882
00883 for (I = Step; I <= Size; I++)
00884 {
00885
00886 Count = -1;
00887 pElement = Matrix->FirstInRow[I];
00888 while (pElement != NULL && pElement->Col < Step)
00889 pElement = pElement->NextInRow;
00890 while (pElement != NULL)
00891 { Count++;
00892 pElement = pElement->NextInRow;
00893 }
00894
00895
00896 ExtRow = Matrix->IntToExtRowMap[I];
00897
00898 #if spSEPARATED_COMPLEX_VECTORS || !spCOMPLEX
00899 if (RHS != NULL)
00900 if (RHS[ExtRow] != 0.0) Count++;
00901 #else
00902 if (RHS != NULL)
00903 { if (Matrix->Complex)
00904 { if ((RHS[2*ExtRow] != 0.0) || (RHS[2*ExtRow+1] != 0.0))
00905 Count++;
00906 }
00907 else if (RHS[I] != 0.0) Count++;
00908 }
00909 #endif
00910 Matrix->MarkowitzRow[I] = Count;
00911 }
00912
00913
00914 for (I = Step; I <= Size; I++)
00915 {
00916
00917 Count = -1;
00918 pElement = Matrix->FirstInCol[I];
00919 while (pElement != NULL && pElement->Row < Step)
00920 pElement = pElement->NextInCol;
00921 while (pElement != NULL)
00922 { Count++;
00923 pElement = pElement->NextInCol;
00924 }
00925 Matrix->MarkowitzCol[I] = Count;
00926 }
00927 return;
00928 }
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967 static void MarkowitzProducts(MatrixPtr Matrix, int Step)
00968 {
00969 int I, *pMarkowitzRow, *pMarkowitzCol;
00970 long Product, *pMarkowitzProduct;
00971 int Size = Matrix->Size;
00972 double fProduct;
00973
00974
00975 Matrix->Singletons = 0;
00976
00977 pMarkowitzProduct = &(Matrix->MarkowitzProd[Step]);
00978 pMarkowitzRow = &(Matrix->MarkowitzRow[Step]);
00979 pMarkowitzCol = &(Matrix->MarkowitzCol[Step]);
00980
00981 for (I = Step; I <= Size; I++)
00982 {
00983
00984 if ((*pMarkowitzRow > SHRT_MAX && *pMarkowitzCol != 0) ||
00985 (*pMarkowitzCol > SHRT_MAX && *pMarkowitzRow != 0))
00986 { fProduct = (double)(*pMarkowitzRow++) * (double)(*pMarkowitzCol++);
00987 if (fProduct >= LONG_MAX)
00988 *pMarkowitzProduct++ = LONG_MAX;
00989 else
00990 *pMarkowitzProduct++ = (long)fProduct;
00991 }
00992 else
00993 { Product = *pMarkowitzRow++ * *pMarkowitzCol++;
00994 if ((*pMarkowitzProduct++ = Product) == 0)
00995 Matrix->Singletons++;
00996 }
00997 }
00998 return;
00999 }
01000
01001 #if DIAGONAL_PIVOTING
01002 #if MODIFIED_MARKOWITZ
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070 static ElementPtr QuicklySearchDiagonal(MatrixPtr Matrix, int Step)
01071 {
01072 long MinMarkowitzProduct, *pMarkowitzProduct;
01073 ElementPtr pDiag, pOtherInRow, pOtherInCol;
01074 int I, NumberOfTies;
01075 ElementPtr ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1];
01076 double Magnitude, LargestInCol, Ratio, MaxRatio;
01077 double LargestOffDiagonal;
01078
01079 NumberOfTies = -1;
01080 MinMarkowitzProduct = LONG_MAX;
01081 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
01082 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
01083
01084
01085 Matrix->MarkowitzProd[Step-1] = -1;
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106 for(;;)
01107 { while (MinMarkowitzProduct < *(--pMarkowitzProduct))
01108 {
01109
01110
01111
01112
01113
01114 }
01115
01116 I = pMarkowitzProduct - Matrix->MarkowitzProd;
01117
01118
01119 if (I < Step) break;
01120 if (I > Matrix->Size) I = Step;
01121
01122 if ((pDiag = Matrix->Diag[I]) == NULL)
01123 continue;
01124 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
01125 continue;
01126
01127 if (*pMarkowitzProduct == 1)
01128 {
01129
01130
01131
01132 pOtherInRow = pDiag->NextInRow;
01133 pOtherInCol = pDiag->NextInCol;
01134 if (pOtherInRow == NULL && pOtherInCol == NULL)
01135 { pOtherInRow = Matrix->FirstInRow[I];
01136 while(pOtherInRow != NULL)
01137 { if (pOtherInRow->Col >= Step && pOtherInRow->Col != I)
01138 break;
01139 pOtherInRow = pOtherInRow->NextInRow;
01140 }
01141 pOtherInCol = Matrix->FirstInCol[I];
01142 while(pOtherInCol != NULL)
01143 { if (pOtherInCol->Row >= Step && pOtherInCol->Row != I)
01144 break;
01145 pOtherInCol = pOtherInCol->NextInCol;
01146 }
01147 }
01148
01149
01150
01151 if (pOtherInRow != NULL && pOtherInCol != NULL)
01152 { if (pOtherInRow->Col == pOtherInCol->Row)
01153 { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
01154 ELEMENT_MAG(pOtherInCol));
01155 if (Magnitude >= LargestOffDiagonal)
01156 {
01157
01158 return pDiag;
01159 }
01160 }
01161 }
01162 }
01163
01164 if (*pMarkowitzProduct < MinMarkowitzProduct)
01165 {
01166
01167 TiedElements[0] = pDiag;
01168 MinMarkowitzProduct = *pMarkowitzProduct;
01169 NumberOfTies = 0;
01170 }
01171 else
01172 {
01173
01174 if (NumberOfTies < MAX_MARKOWITZ_TIES)
01175 { TiedElements[++NumberOfTies] = pDiag;
01176 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
01177 break;
01178 }
01179 }
01180 }
01181
01182
01183 if (NumberOfTies < 0)
01184 return NULL;
01185
01186
01187 ChosenPivot = NULL;
01188 MaxRatio = 1.0 / Matrix->RelThreshold;
01189
01190 for (I = 0; I <= NumberOfTies; I++)
01191 { pDiag = TiedElements[I];
01192 Magnitude = ELEMENT_MAG(pDiag);
01193 LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
01194 Ratio = LargestInCol / Magnitude;
01195 if (Ratio < MaxRatio)
01196 { ChosenPivot = pDiag;
01197 MaxRatio = Ratio;
01198 }
01199 }
01200 return ChosenPivot;
01201 }
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212 #else
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263 static ElementPtr QuicklySearchDiagonal(MatrixPtr Matrix, int Step)
01264 {
01265 long MinMarkowitzProduct, *pMarkowitzProduct;
01266 ElementPtr pDiag;
01267 int I;
01268 ElementPtr ChosenPivot, pOtherInRow, pOtherInCol;
01269 double Magnitude, LargestInCol, LargestOffDiagonal;
01270
01271
01272 ChosenPivot = NULL;
01273 MinMarkowitzProduct = LONG_MAX;
01274 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
01275 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
01276
01277
01278 Matrix->MarkowitzProd[Step-1] = -1;
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299 for (;;)
01300 { while (*(--pMarkowitzProduct) >= MinMarkowitzProduct)
01301 {
01302 }
01303
01304 I = (int)(pMarkowitzProduct - Matrix->MarkowitzProd);
01305
01306
01307 if (I < Step) break;
01308 if (I > Matrix->Size) I = Step;
01309
01310 if ((pDiag = Matrix->Diag[I]) == NULL)
01311 continue;
01312 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
01313 continue;
01314
01315 if (*pMarkowitzProduct == 1)
01316 {
01317
01318
01319
01320 pOtherInRow = pDiag->NextInRow;
01321 pOtherInCol = pDiag->NextInCol;
01322 if (pOtherInRow == NULL && pOtherInCol == NULL)
01323 { pOtherInRow = Matrix->FirstInRow[I];
01324 while(pOtherInRow != NULL)
01325 { if (pOtherInRow->Col >= Step && pOtherInRow->Col != I)
01326 break;
01327 pOtherInRow = pOtherInRow->NextInRow;
01328 }
01329 pOtherInCol = Matrix->FirstInCol[I];
01330 while(pOtherInCol != NULL)
01331 { if (pOtherInCol->Row >= Step && pOtherInCol->Row != I)
01332 break;
01333 pOtherInCol = pOtherInCol->NextInCol;
01334 }
01335 }
01336
01337
01338
01339 if (pOtherInRow != NULL && pOtherInCol != NULL)
01340 { if (pOtherInRow->Col == pOtherInCol->Row)
01341 { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
01342 ELEMENT_MAG(pOtherInCol));
01343 if (Magnitude >= LargestOffDiagonal)
01344 {
01345
01346 return pDiag;
01347 }
01348 }
01349 }
01350 }
01351
01352 MinMarkowitzProduct = *pMarkowitzProduct;
01353 ChosenPivot = pDiag;
01354 }
01355
01356 if (ChosenPivot != NULL)
01357 { LargestInCol = FindBiggestInColExclude( Matrix, ChosenPivot, Step );
01358 if( ELEMENT_MAG(ChosenPivot) <= Matrix->RelThreshold * LargestInCol )
01359 ChosenPivot = NULL;
01360 }
01361 return ChosenPivot;
01362 }
01363 #endif
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423 static ElementPtr SearchDiagonal(MatrixPtr Matrix, int Step)
01424 {
01425 int J;
01426 long MinMarkowitzProduct, *pMarkowitzProduct;
01427 int I;
01428 ElementPtr pDiag;
01429 int NumberOfTies, Size = Matrix->Size;
01430 ElementPtr ChosenPivot;
01431 double Magnitude, Ratio, RatioOfAccepted, LargestInCol;
01432
01433
01434 ChosenPivot = NULL;
01435 MinMarkowitzProduct = LONG_MAX;
01436 pMarkowitzProduct = &(Matrix->MarkowitzProd[Size+2]);
01437 Matrix->MarkowitzProd[Size+1] = Matrix->MarkowitzProd[Step];
01438
01439
01440 for (J = Size+1; J > Step; J--)
01441 {
01442 if (*(--pMarkowitzProduct) > MinMarkowitzProduct)
01443 continue;
01444 if (J > Matrix->Size)
01445 I = Step;
01446 else
01447 I = J;
01448 if ((pDiag = Matrix->Diag[I]) == NULL)
01449 continue;
01450 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
01451 continue;
01452
01453
01454 LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
01455 if (Magnitude <= Matrix->RelThreshold * LargestInCol)
01456 continue;
01457
01458 if (*pMarkowitzProduct < MinMarkowitzProduct)
01459 {
01460
01461 ChosenPivot = pDiag;
01462 MinMarkowitzProduct = *pMarkowitzProduct;
01463 RatioOfAccepted = LargestInCol / Magnitude;
01464 NumberOfTies = 0;
01465 }
01466 else
01467 {
01468
01469 NumberOfTies++;
01470 Ratio = LargestInCol / Magnitude;
01471 if (Ratio < RatioOfAccepted)
01472 { ChosenPivot = pDiag;
01473 RatioOfAccepted = Ratio;
01474 }
01475 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
01476 return ChosenPivot;
01477 }
01478 }
01479 return ChosenPivot;
01480 }
01481 #endif
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538 static ElementPtr SearchEntireMatrix(MatrixPtr Matrix, int Step)
01539 {
01540 int I, Size = Matrix->Size;
01541 ElementPtr pElement;
01542 int NumberOfTies;
01543 long Product, MinMarkowitzProduct;
01544 ElementPtr ChosenPivot, pLargestElement;
01545 double Magnitude, LargestElementMag, Ratio, RatioOfAccepted, LargestInCol;
01546
01547
01548 ChosenPivot = NULL;
01549 LargestElementMag = 0.0;
01550 MinMarkowitzProduct = LONG_MAX;
01551
01552
01553 for (I = Step; I <= Size; I++)
01554 { pElement = Matrix->FirstInCol[I];
01555
01556 while (pElement != NULL && pElement->Row < Step)
01557 pElement = pElement->NextInCol;
01558
01559 if((LargestInCol = FindLargestInCol(pElement)) == 0.0)
01560 continue;
01561
01562 while (pElement != NULL)
01563 {
01564
01565
01566 if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElementMag)
01567 { LargestElementMag = Magnitude;
01568 pLargestElement = pElement;
01569 }
01570
01571 Product = Matrix->MarkowitzRow[pElement->Row] *
01572 Matrix->MarkowitzCol[pElement->Col];
01573
01574
01575 if ((Product <= MinMarkowitzProduct) &&
01576 (Magnitude > Matrix->RelThreshold * LargestInCol) &&
01577 (Magnitude > Matrix->AbsThreshold))
01578 {
01579
01580
01581 if (Product < MinMarkowitzProduct)
01582 {
01583
01584 ChosenPivot = pElement;
01585 MinMarkowitzProduct = Product;
01586 RatioOfAccepted = LargestInCol / Magnitude;
01587 NumberOfTies = 0;
01588 }
01589 else
01590 {
01591
01592 NumberOfTies++;
01593 Ratio = LargestInCol / Magnitude;
01594 if (Ratio < RatioOfAccepted)
01595 { ChosenPivot = pElement;
01596 RatioOfAccepted = Ratio;
01597 }
01598 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
01599 return ChosenPivot;
01600 }
01601 }
01602 pElement = pElement->NextInCol;
01603 }
01604 }
01605
01606 if (ChosenPivot != NULL) return ChosenPivot;
01607
01608 if (LargestElementMag == 0.0)
01609 { Matrix->Error = spSINGULAR;
01610 return NULL;
01611 }
01612
01613 Matrix->Error = spSMALL_PIVOT;
01614 return pLargestElement;
01615 }
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656 static double FindLargestInCol(ElementPtr pElement)
01657 {
01658 double Magnitude, Largest = 0.0;
01659
01660
01661
01662 while (pElement != NULL)
01663 { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
01664 Largest = Magnitude;
01665 pElement = pElement->NextInCol;
01666 }
01667
01668 return Largest;
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
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719 static double FindBiggestInColExclude(MatrixPtr Matrix, ElementPtr pElement, int Step)
01720 {
01721 int Row, Col;
01722 double Largest, Magnitude;
01723
01724
01725 Row = pElement->Row;
01726 Col = pElement->Col;
01727 pElement = Matrix->FirstInCol[Col];
01728
01729
01730 while ((pElement != NULL) && (pElement->Row < Step))
01731 pElement = pElement->NextInCol;
01732
01733
01734 if (pElement->Row != Row)
01735 Largest = ELEMENT_MAG(pElement);
01736 else
01737 Largest = 0.0;
01738
01739
01740 while ((pElement = pElement->NextInCol) != NULL)
01741 { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
01742 { if (pElement->Row != Row)
01743 Largest = Magnitude;
01744 }
01745 }
01746
01747 return Largest;
01748 }
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789 void spcRowExchange(MatrixPtr Matrix, int Row1, int Row2 );
01790 static void ExchangeColElements(MatrixPtr Matrix, int Row1, ElementPtr Element1, int Row2, ElementPtr Element2, int Column );
01791 static void ExchangeRowElements(MatrixPtr Matrix, int Col1, ElementPtr Element1, int Col2, ElementPtr Element2, int Row );
01792 ElementPtr spcFindElementInCol(MatrixPtr Matrix, ElementPtr *LastAddr, int Row, int Col, BOOLEAN CreateIfMissing);
01793
01794 static void ExchangeRowsAndCols(MatrixPtr Matrix, ElementPtr pPivot, int Step )
01795 {
01796 int Row, Col;
01797 long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col;
01798
01799 Row = pPivot->Row;
01800 Col = pPivot->Col;
01801 Matrix->PivotsOriginalRow = Row;
01802 Matrix->PivotsOriginalCol = Col;
01803
01804 if ((Row == Step) && (Col == Step)) return;
01805
01806
01807 if (Row == Col)
01808 { spcRowExchange( Matrix, Step, Row );
01809 spcColExchange( Matrix, Step, Col );
01810 SWAP( long, Matrix->MarkowitzProd[Step], Matrix->MarkowitzProd[Row] );
01811 SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Step] );
01812 }
01813 else
01814 {
01815
01816
01817 OldMarkowitzProd_Step = Matrix->MarkowitzProd[Step];
01818 OldMarkowitzProd_Row = Matrix->MarkowitzProd[Row];
01819 OldMarkowitzProd_Col = Matrix->MarkowitzProd[Col];
01820
01821
01822 if (Row != Step)
01823 { spcRowExchange( Matrix, Step, Row );
01824 Matrix->NumberOfInterchangesIsOdd = !Matrix->NumberOfInterchangesIsOdd;
01825 Matrix->MarkowitzProd[Row] = Matrix->MarkowitzRow[Row] *
01826 Matrix->MarkowitzCol[Row];
01827
01828
01829 if ((Matrix->MarkowitzProd[Row]==0) != (OldMarkowitzProd_Row==0))
01830 { if (OldMarkowitzProd_Row == 0)
01831 Matrix->Singletons--;
01832 else
01833 Matrix->Singletons++;
01834 }
01835 }
01836
01837
01838 if (Col != Step)
01839 { spcColExchange( Matrix, Step, Col );
01840 Matrix->NumberOfInterchangesIsOdd =
01841 !Matrix->NumberOfInterchangesIsOdd;
01842 Matrix->MarkowitzProd[Col] = Matrix->MarkowitzCol[Col] *
01843 Matrix->MarkowitzRow[Col];
01844
01845
01846 if ((Matrix->MarkowitzProd[Col]==0) != (OldMarkowitzProd_Col==0))
01847 { if (OldMarkowitzProd_Col == 0)
01848 Matrix->Singletons--;
01849 else
01850 Matrix->Singletons++;
01851 }
01852
01853 Matrix->Diag[Col] = spcFindElementInCol( Matrix,
01854 Matrix->FirstInCol+Col,
01855 Col, Col, 0 );
01856 }
01857 if (Row != Step)
01858 { Matrix->Diag[Row] = spcFindElementInCol( Matrix,
01859 Matrix->FirstInCol+Row,
01860 Row, Row, 0 );
01861 }
01862 Matrix->Diag[Step] = spcFindElementInCol( Matrix,
01863 Matrix->FirstInCol+Step,
01864 Step, Step, 0 );
01865
01866
01867 Matrix->MarkowitzProd[Step] = Matrix->MarkowitzCol[Step] *
01868 Matrix->MarkowitzRow[Step];
01869 if ((Matrix->MarkowitzProd[Step]==0) != (OldMarkowitzProd_Step==0))
01870 { if (OldMarkowitzProd_Step == 0)
01871 Matrix->Singletons--;
01872 else
01873 Matrix->Singletons++;
01874 }
01875 }
01876 return;
01877 }
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916 void spcRowExchange(MatrixPtr Matrix, int Row1, int Row2 )
01917 {
01918 ElementPtr Row1Ptr, Row2Ptr, Element1, Element2;
01919 int Column;
01920
01921 if (Row1 > Row2) SWAP(int, Row1, Row2);
01922
01923 Row1Ptr = Matrix->FirstInRow[Row1];
01924 Row2Ptr = Matrix->FirstInRow[Row2];
01925 while (Row1Ptr != NULL || Row2Ptr != NULL)
01926 {
01927
01928 if (Row1Ptr == NULL)
01929 {
01930 Column = Row2Ptr->Col;
01931 Element1 = NULL;
01932 Element2 = Row2Ptr;
01933 Row2Ptr = Row2Ptr->NextInRow;
01934 }
01935 else if (Row2Ptr == NULL)
01936 {
01937 Column = Row1Ptr->Col;
01938 Element1 = Row1Ptr;
01939 Element2 = NULL;
01940 Row1Ptr = Row1Ptr->NextInRow;
01941 }
01942 else if (Row1Ptr->Col < Row2Ptr->Col)
01943 {
01944 Column = Row1Ptr->Col;
01945 Element1 = Row1Ptr;
01946 Element2 = NULL;
01947 Row1Ptr = Row1Ptr->NextInRow;
01948 }
01949 else if (Row1Ptr->Col > Row2Ptr->Col)
01950 {
01951 Column = Row2Ptr->Col;
01952 Element1 = NULL;
01953 Element2 = Row2Ptr;
01954 Row2Ptr = Row2Ptr->NextInRow;
01955 }
01956 else
01957 {
01958 Column = Row1Ptr->Col;
01959 Element1 = Row1Ptr;
01960 Element2 = Row2Ptr;
01961 Row1Ptr = Row1Ptr->NextInRow;
01962 Row2Ptr = Row2Ptr->NextInRow;
01963 }
01964 ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column);
01965 }
01966
01967 if (Matrix->InternalVectorsAllocated)
01968 SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]);
01969 SWAP( ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]);
01970 SWAP( int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]);
01971 #if TRANSLATE
01972 Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row1] ] = Row1;
01973 Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row2] ] = Row2;
01974 #endif
01975
01976 return;
01977 }
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006 void spcColExchange(MatrixPtr Matrix, int Col1, int Col2 )
02007 {
02008 ElementPtr Col1Ptr, Col2Ptr, Element1, Element2;
02009 int Row;
02010
02011 if (Col1 > Col2) SWAP(int, Col1, Col2);
02012
02013 Col1Ptr = Matrix->FirstInCol[Col1];
02014 Col2Ptr = Matrix->FirstInCol[Col2];
02015 while (Col1Ptr != NULL || Col2Ptr != NULL)
02016 {
02017
02018 if (Col1Ptr == NULL)
02019 { Row = Col2Ptr->Row;
02020 Element1 = NULL;
02021 Element2 = Col2Ptr;
02022 Col2Ptr = Col2Ptr->NextInCol;
02023 }
02024 else if (Col2Ptr == NULL)
02025 { Row = Col1Ptr->Row;
02026 Element1 = Col1Ptr;
02027 Element2 = NULL;
02028 Col1Ptr = Col1Ptr->NextInCol;
02029 }
02030 else if (Col1Ptr->Row < Col2Ptr->Row)
02031 { Row = Col1Ptr->Row;
02032 Element1 = Col1Ptr;
02033 Element2 = NULL;
02034 Col1Ptr = Col1Ptr->NextInCol;
02035 }
02036 else if (Col1Ptr->Row > Col2Ptr->Row)
02037 { Row = Col2Ptr->Row;
02038 Element1 = NULL;
02039 Element2 = Col2Ptr;
02040 Col2Ptr = Col2Ptr->NextInCol;
02041 }
02042 else
02043 { Row = Col1Ptr->Row;
02044 Element1 = Col1Ptr;
02045 Element2 = Col2Ptr;
02046 Col1Ptr = Col1Ptr->NextInCol;
02047 Col2Ptr = Col2Ptr->NextInCol;
02048 }
02049
02050 ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row);
02051 }
02052
02053 if (Matrix->InternalVectorsAllocated)
02054 SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]);
02055 SWAP( ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
02056 SWAP( int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
02057 #if TRANSLATE
02058 Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col1] ] = Col1;
02059 Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col2] ] = Col2;
02060 #endif
02061 return;
02062 }
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099 static void ExchangeColElements(MatrixPtr Matrix, int Row1, ElementPtr Element1, int Row2, ElementPtr Element2, int Column )
02100 {
02101 ElementPtr *ElementAboveRow1, *ElementAboveRow2, ElementBelowRow1, ElementBelowRow2, pElement;
02102
02103
02104 ElementAboveRow1 = &(Matrix->FirstInCol[Column]);
02105 pElement = *ElementAboveRow1;
02106 while (pElement->Row < Row1)
02107 {
02108 ElementAboveRow1 = &(pElement->NextInCol);
02109 pElement = *ElementAboveRow1;
02110 }
02111 if (Element1 != NULL)
02112 {
02113 ElementBelowRow1 = Element1->NextInCol;
02114 if (Element2 == NULL)
02115 {
02116
02117 if ( ElementBelowRow1 != NULL && ElementBelowRow1->Row < Row2 )
02118 {
02119
02120 *ElementAboveRow1 = ElementBelowRow1;
02121
02122
02123 pElement = ElementBelowRow1;
02124 do
02125 {
02126 ElementAboveRow2 = &(pElement->NextInCol);
02127 pElement = *ElementAboveRow2;
02128 } while (pElement != NULL && pElement->Row < Row2);
02129
02130
02131 *ElementAboveRow2 = Element1;
02132 Element1->NextInCol = pElement;
02133 *ElementAboveRow1 =ElementBelowRow1;
02134 }
02135 Element1->Row = Row2;
02136 }
02137 else
02138 {
02139
02140 if ( ElementBelowRow1->Row == Row2)
02141 {
02142
02143 Element1->NextInCol = Element2->NextInCol;
02144 Element2->NextInCol = Element1;
02145 *ElementAboveRow1 = Element2;
02146 }
02147 else
02148 {
02149
02150 pElement = ElementBelowRow1;
02151 do
02152 {
02153 ElementAboveRow2 = &(pElement->NextInCol);
02154 pElement = *ElementAboveRow2;
02155 } while (pElement->Row < Row2);
02156
02157 ElementBelowRow2 = Element2->NextInCol;
02158
02159
02160 *ElementAboveRow1 = Element2;
02161 Element2->NextInCol = ElementBelowRow1;
02162 *ElementAboveRow2 = Element1;
02163 Element1->NextInCol = ElementBelowRow2;
02164 }
02165 Element1->Row = Row2;
02166 Element2->Row = Row1;
02167 }
02168 }
02169 else
02170 {
02171
02172 ElementBelowRow1 = pElement;
02173
02174
02175 if (ElementBelowRow1->Row != Row2)
02176 {
02177 do
02178 {
02179 ElementAboveRow2 = &(pElement->NextInCol);
02180 pElement = *ElementAboveRow2;
02181 } while (pElement->Row < Row2);
02182 ElementBelowRow2 = Element2->NextInCol;
02183
02184
02185 *ElementAboveRow2 = Element2->NextInCol;
02186 *ElementAboveRow1 = Element2;
02187 Element2->NextInCol = ElementBelowRow1;
02188 }
02189 Element2->Row = Row1;
02190 }
02191 return;
02192 }
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231 static void ExchangeRowElements(MatrixPtr Matrix, int Col1, ElementPtr Element1, int Col2, ElementPtr Element2, int Row )
02232 {
02233 ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2, ElementRightOfCol1, ElementRightOfCol2, pElement;
02234
02235
02236 ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]);
02237 pElement = *ElementLeftOfCol1;
02238 while (pElement->Col < Col1)
02239 {
02240 ElementLeftOfCol1 = &(pElement->NextInRow);
02241 pElement = *ElementLeftOfCol1;
02242 }
02243 if (Element1 != NULL)
02244 {
02245 ElementRightOfCol1 = Element1->NextInRow;
02246 if (Element2 == NULL)
02247 {
02248
02249 if ( ElementRightOfCol1 != NULL && ElementRightOfCol1->Col < Col2 )
02250 {
02251
02252 *ElementLeftOfCol1 = ElementRightOfCol1;
02253
02254
02255 pElement = ElementRightOfCol1;
02256 do
02257 {
02258 ElementLeftOfCol2 = &(pElement->NextInRow);
02259 pElement = *ElementLeftOfCol2;
02260 } while (pElement != NULL && pElement->Col < Col2);
02261
02262
02263 *ElementLeftOfCol2 = Element1;
02264 Element1->NextInRow = pElement;
02265 *ElementLeftOfCol1 =ElementRightOfCol1;
02266 }
02267 Element1->Col = Col2;
02268 }
02269 else
02270 {
02271
02272 if ( ElementRightOfCol1->Col == Col2)
02273 {
02274
02275 Element1->NextInRow = Element2->NextInRow;
02276 Element2->NextInRow = Element1;
02277 *ElementLeftOfCol1 = Element2;
02278 }
02279 else
02280 {
02281
02282 pElement = ElementRightOfCol1;
02283 do
02284 {
02285 ElementLeftOfCol2 = &(pElement->NextInRow);
02286 pElement = *ElementLeftOfCol2;
02287 } while (pElement->Col < Col2);
02288
02289 ElementRightOfCol2 = Element2->NextInRow;
02290
02291
02292 *ElementLeftOfCol1 = Element2;
02293 Element2->NextInRow = ElementRightOfCol1;
02294 *ElementLeftOfCol2 = Element1;
02295 Element1->NextInRow = ElementRightOfCol2;
02296 }
02297 Element1->Col = Col2;
02298 Element2->Col = Col1;
02299 }
02300 }
02301 else
02302 {
02303
02304 ElementRightOfCol1 = pElement;
02305
02306
02307 if (ElementRightOfCol1->Col != Col2)
02308 {
02309 do
02310 {
02311 ElementLeftOfCol2 = &(pElement->NextInRow);
02312 pElement = *ElementLeftOfCol2;
02313 } while (pElement->Col < Col2);
02314
02315 ElementRightOfCol2 = Element2->NextInRow;
02316
02317
02318 *ElementLeftOfCol2 = Element2->NextInRow;
02319 *ElementLeftOfCol1 = Element2;
02320 Element2->NextInRow = ElementRightOfCol1;
02321 }
02322 Element2->Col = Col1;
02323 }
02324 return;
02325 }
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353 static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot)
02354 {
02355 #if REAL
02356 ElementPtr pSub, pLower, pUpper;
02357 int Row;
02358
02359
02360
02361
02362 if (ABSS(pPivot->Real) == 0.0)
02363 { (void)MatrixIsSingular( Matrix, pPivot->Row );
02364 return;
02365 }
02366 pPivot->Real = 1.0 / pPivot->Real;
02367
02368 pUpper = pPivot->NextInRow;
02369 while (pUpper != NULL)
02370 {
02371
02372 pUpper->Real *= pPivot->Real;
02373
02374 pSub = pUpper->NextInCol;
02375 pLower = pPivot->NextInCol;
02376 while (pLower != NULL)
02377 { Row = pLower->Row;
02378
02379
02380 while (pSub != NULL && pSub->Row < Row)
02381 pSub = pSub->NextInCol;
02382
02383
02384 if (pSub == NULL || pSub->Row > Row)
02385 { pSub = CreateFillin( Matrix, Row, pUpper->Col );
02386 if (pSub == NULL)
02387 { Matrix->Error = spNO_MEMORY;
02388 return;
02389 }
02390 }
02391 pSub->Real -= pUpper->Real * pLower->Real;
02392 pSub = pSub->NextInCol;
02393 pLower = pLower->NextInCol;
02394 }
02395 pUpper = pUpper->NextInRow;
02396 }
02397 return;
02398 #endif
02399 }
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423 static void UpdateMarkowitzNumbers(MatrixPtr Matrix, ElementPtr pPivot)
02424 {
02425 int Row, Col;
02426 ElementPtr ColPtr, RowPtr;
02427 int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
02428 double Product;
02429
02430
02431
02432
02433 for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol)
02434 { Row = ColPtr->Row;
02435 --MarkoRow[Row];
02436
02437
02438 if ((MarkoRow[Row] > SHRT_MAX && MarkoCol[Row] != 0) ||
02439 (MarkoCol[Row] > SHRT_MAX && MarkoRow[Row] != 0))
02440 { Product = MarkoCol[Row] * MarkoRow[Row];
02441 if (Product >= LONG_MAX)
02442 Matrix->MarkowitzProd[Row] = LONG_MAX;
02443 else
02444 Matrix->MarkowitzProd[Row] = (long)Product;
02445 }
02446 else Matrix->MarkowitzProd[Row] = MarkoRow[Row] * MarkoCol[Row];
02447 if (MarkoRow[Row] == 0)
02448 Matrix->Singletons++;
02449 }
02450
02451 for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow)
02452 { Col = RowPtr->Col;
02453 --MarkoCol[Col];
02454
02455
02456 if ((MarkoRow[Col] > SHRT_MAX && MarkoCol[Col] != 0) ||
02457 (MarkoCol[Col] > SHRT_MAX && MarkoRow[Col] != 0))
02458 { Product = MarkoCol[Col] * MarkoRow[Col];
02459 if (Product >= LONG_MAX)
02460 Matrix->MarkowitzProd[Col] = LONG_MAX;
02461 else
02462 Matrix->MarkowitzProd[Col] = (long)Product;
02463 }
02464 else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col];
02465 if ((MarkoCol[Col] == 0) && (MarkoRow[Col] != 0))
02466 Matrix->Singletons++;
02467 }
02468 return;
02469 }
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500 static ElementPtr CreateFillin(MatrixPtr Matrix, int Row, int Col)
02501 {
02502 ElementPtr pElement, *ppElementAbove;
02503
02504
02505
02506
02507 ppElementAbove = &Matrix->FirstInCol[Col];
02508 pElement = *ppElementAbove;
02509 while (pElement != NULL)
02510 { if (pElement->Row < Row)
02511 { ppElementAbove = &pElement->NextInCol;
02512 pElement = *ppElementAbove;
02513 }
02514 else break;
02515 }
02516
02517
02518 pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, 1);
02519
02520
02521 Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] *
02522 Matrix->MarkowitzCol[Row];
02523 if ((Matrix->MarkowitzRow[Row] == 1) && (Matrix->MarkowitzCol[Row] != 0))
02524 Matrix->Singletons--;
02525 Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] *
02526 Matrix->MarkowitzRow[Col];
02527 if ((Matrix->MarkowitzRow[Col] != 0) && (Matrix->MarkowitzCol[Col] == 1))
02528 Matrix->Singletons--;
02529
02530 return pElement;
02531 }
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547 static int MatrixIsSingular(MatrixPtr Matrix, int Step)
02548 {
02549 Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
02550 Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
02551 return (Matrix->Error = spSINGULAR);
02552 }
02553
02554 static int ZeroPivot(MatrixPtr Matrix, int Step)
02555 {
02556 Matrix->SingularRow = Matrix->IntToExtRowMap[Step];
02557 Matrix->SingularCol = Matrix->IntToExtColMap[Step];
02558 return Matrix->Error = spZERO_DIAG;
02559 }
02560 #if ANNOTATE == FULL
02561
02562
02563
02564
02565
02566
02567
02568 static WriteStatus(MatrixPtr Matrix, int Step)
02569 {
02570 int I;
02571
02572 printf("Step = %1d ", Step);
02573 printf("Pivot found at %1d,%1d using ", Matrix->PivotsOriginalRow,
02574 Matrix->PivotsOriginalCol);
02575 switch(Matrix->PivotSelectionMethod)
02576 {
02577 case 's': printf("SearchForSingleton\n"); break;
02578 case 'q': printf("QuicklySearchDiagonal\n"); break;
02579 case 'd': printf("SearchDiagonal\n"); break;
02580 case 'e': printf("SearchEntireMatrix\n"); break;
02581 }
02582
02583 printf("MarkowitzRow = ");
02584 for (I = 1; I <= Matrix->Size; I++)
02585 printf("%2d ", Matrix->MarkowitzRow[I]);
02586 printf("\n");
02587
02588 printf("MarkowitzCol = ");
02589 for (I = 1; I <= Matrix->Size; I++)
02590 printf("%2d ", Matrix->MarkowitzCol[I]);
02591 printf("\n");
02592
02593 printf("MarkowitzProduct = ");
02594 for (I = 1; I <= Matrix->Size; I++)
02595 printf("%2d ", Matrix->MarkowitzProd[I]);
02596 printf("\n");
02597
02598 printf("Singletons = %2d\n", Matrix->Singletons);
02599
02600 printf("IntToExtRowMap = ");
02601 for (I = 1; I <= Matrix->Size; I++)
02602 printf("%2d ", Matrix->IntToExtRowMap[I]);
02603 printf("\n");
02604
02605 printf("IntToExtColMap = ");
02606 for (I = 1; I <= Matrix->Size; I++)
02607 printf("%2d ", Matrix->IntToExtColMap[I]);
02608 printf("\n");
02609
02610 printf("ExtToIntRowMap = ");
02611 for (I = 1; I <= Matrix->ExtSize; I++)
02612 printf("%2d ", Matrix->ExtToIntRowMap[I]);
02613 printf("\n");
02614
02615 printf("ExtToIntColMap = ");
02616 for (I = 1; I <= Matrix->ExtSize; I++)
02617 printf("%2d ", Matrix->ExtToIntColMap[I]);
02618 printf("\n\n");
02619
02620 return;
02621 }
02622 #endif