1 /*
2 * MATRIX FACTORIZATION MODULE
3 *
4 * Author: Advising Professor:
5 * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
6 * UC Berkeley
7 *
8 * This file contains the routines to factor the matrix into LU form.
9 *
10 * >>> User accessible functions contained in this file:
11 * spOrderAndFactor
12 * spFactor
13 * spPartition
14 *
15 * >>> Other functions contained in this file:
16 * spcRowExchange
17 * spcColExchange
18 */
19
20
21 /*
22 * Revision and copyright information.
23 *
24 * Copyright (c) 1985,86,87,88
25 * by Kenneth S. Kundert and the University of California.
26 *
27 * Permission to use, copy, modify, and distribute this software and
28 * its documentation for any purpose and without fee is hereby granted,
29 * provided that the copyright notices appear in all copies and
30 * supporting documentation and that the authors and the University of
31 * California are properly credited. The authors and the University of
32 * California make no representations as to the suitability of this
33 * software for any purpose. It is provided `as is', without express
34 * or implied warranty.
35 */
36
37 #ifndef lint
38 static char copyright[] =
39 "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
40 static char RCSid[] =
41 "@(#)$Header: spFactor.c,v 1.3 88/06/24 05:01:12 kundert Exp $";
42 #endif
43
44
45
46 /*
47 * IMPORTS
48 *
49 * >>> Import descriptions:
50 * spConfig.h
51 * Macros that customize the sparse matrix routines.
52 * spMatrix.h
53 * Macros and declarations to be imported by the user.
54 * spDefs.h
55 * Matrix type and macro definitions for the sparse matrix routines.
56 */
57
58 #define spINSIDE_SPARSE
59 #include "spConfig.h"
60 #include "spMatrix.h"
61 #include "spDefs.h"
62
63 static FactorComplexMatrix();
64 static CreateInternalVectors();
65 static CountMarkowitz();
66 static MarkowitzProducts();
67 static ElementPtr SearchForPivot();
68 static ElementPtr SearchForSingleton();
69 static ElementPtr QuicklySearchDiagonal();
70 static ElementPtr SearchDiagonal();
71 static ElementPtr SearchEntireMatrix();
72 static RealNumber FindLargestInCol();
73 static RealNumber FindBiggestInColExclude();
74 static ExchangeRowsAndCols();
75 static ExchangeColElements();
76 static ExchangeRowElements();
77 static RealRowColElimination();
78 static ComplexRowColElimination();
79 static UpdateMarkowitzNumbers();
80 static ElementPtr CreateFillin();
81 static MatrixIsSingular();
82 static ZeroPivot();
83 static WriteStatus();
84
85
86
87
88
89 /*
90 * ORDER AND FACTOR MATRIX
91 *
92 * This routine chooses a pivot order for the matrix and factors it
93 * into LU form. It handles both the initial factorization and subsequent
94 * factorizations when a reordering is desired. This is handled in a manner
95 * that is transparent to the user. The routine uses a variation of
96 * Gauss's method where the pivots are associated with L and the
97 * diagonal terms of U are one.
98 *
99 * >>> Returned:
100 * The error code is returned. Possible errors are listed below.
101 *
102 * >>> Arguments:
103 * Matrix <input> (char *)
104 * Pointer to matrix.
105 * RHS <input> (RealVector)
106 * Representative right-hand side vector that is used to determine
107 * pivoting order when the right hand side vector is sparse. If
108 * RHS is a NULL pointer then the RHS vector is assumed to
109 * be full and it is not used when determining the pivoting
110 * order.
111 * RelThreshold <input> (RealNumber)
112 * This number determines what the pivot relative threshold will
113 * be. It should be between zero and one. If it is one then the
114 * pivoting method becomes complete pivoting, which is very slow
115 * and tends to fill up the matrix. If it is set close to zero
116 * the pivoting method becomes strict Markowitz with no
117 * threshold. The pivot threshold is used to eliminate pivot
118 * candidates that would cause excessive element growth if they
119 * were used. Element growth is the cause of roundoff error.
120 * Element growth occurs even in well-conditioned matrices.
121 * Setting the RelThreshold large will reduce element growth and
122 * roundoff error, but setting it too large will cause execution
123 * time to be excessive and will result in a large number of
124 * fill-ins. If this occurs, accuracy can actually be degraded
125 * because of the large number of operations required on the
126 * matrix due to the large number of fill-ins. A good value seems
127 * to be 0.001. The default is chosen by giving a value larger
128 * than one or less than or equal to zero. This value should be
129 * increased and the matrix resolved if growth is found to be
130 * excessive. Changing the pivot threshold does not improve
131 * performance on matrices where growth is low, as is often the
132 * case with ill-conditioned matrices. Once a valid threshold is
133 * given, it becomes the new default. The default value of
134 * RelThreshold was choosen for use with nearly diagonally
135 * dominant matrices such as node- and modified-node admittance
136 * matrices. For these matrices it is usually best to use
137 * diagonal pivoting. For matrices without a strong diagonal, it
138 * is usually best to use a larger threshold, such as 0.01 or
139 * 0.1.
140 * AbsThreshold <input> (RealNumber)
141 * The absolute magnitude an element must have to be considered
142 * as a pivot candidate, except as a last resort. This number
143 * should be set significantly smaller than the smallest diagonal
144 * element that is is expected to be placed in the matrix. If
145 * there is no reasonable prediction for the lower bound on these
146 * elements, then AbsThreshold should be set to zero.
147 * AbsThreshold is used to reduce the possibility of choosing as a
148 * pivot an element that has suffered heavy cancellation and as a
149 * result mainly consists of roundoff error. Once a valid
150 * threshold is given, it becomes the new default.
151 * DiagPivoting <input> (BOOLEAN)
152 * A flag indicating that pivot selection should be confined to the
153 * diagonal if possible. If DiagPivoting is nonzero and if
154 * DIAGONAL_PIVOTING is enabled pivots will be chosen only from
155 * the diagonal unless there are no diagonal elements that satisfy
156 * the threshold criteria. Otherwise, the entire reduced
157 * submatrix is searched when looking for a pivot. The diagonal
158 * pivoting in Sparse is efficient and well refined, while the
159 * off-diagonal pivoting is not. For symmetric and near symmetric
160 * matrices, it is best to use diagonal pivoting because it
161 * results in the best performance when reordering the matrix and
162 * when factoring the matrix without ordering. If there is a
163 * considerable amount of nonsymmetry in the matrix, then
164 * off-diagonal pivoting may result in a better equation ordering
165 * simply because there are more pivot candidates to choose from.
166 * A better ordering results in faster subsequent factorizations.
167 * However, the initial pivot selection process takes considerably
168 * longer for off-diagonal pivoting.
169 *
170 * >>> Local variables:
171 * pPivot (ElementPtr)
172 * Pointer to the element being used as a pivot.
173 * ReorderingRequired (BOOLEAN)
174 * Flag that indicates whether reordering is required.
175 *
176 * >>> Possible errors:
177 * spNO_MEMORY
178 * spSINGULAR
179 * spSMALL_PIVOT
180 * Error is cleared in this function.
181 */
182
183 int
spOrderAndFactor(eMatrix,RHS,RelThreshold,AbsThreshold,DiagPivoting)184 spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )
185
186 char *eMatrix;
187 RealNumber RHS[], RelThreshold, AbsThreshold;
188 BOOLEAN DiagPivoting;
189 {
190 MatrixPtr Matrix = (MatrixPtr)eMatrix;
191 ElementPtr pPivot;
192 int Step, Size, ReorderingRequired;
193 ElementPtr SearchForPivot();
194 RealNumber LargestInCol, FindLargestInCol();
195
196 /* Begin `spOrderAndFactor'. */
197 ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
198
199 Matrix->Error = spOKAY;
200 Size = Matrix->Size;
201 if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold;
202 if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold;
203 Matrix->RelThreshold = RelThreshold;
204 if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold;
205 Matrix->AbsThreshold = AbsThreshold;
206 ReorderingRequired = NO;
207
208 if (NOT Matrix->NeedsOrdering)
209 {
210 /* Matrix has been factored before and reordering is not required. */
211 for (Step = 1; Step <= Size; Step++)
212 { pPivot = Matrix->Diag[Step];
213 LargestInCol = FindLargestInCol(pPivot->NextInCol);
214 if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
215 { if (Matrix->Complex)
216 ComplexRowColElimination( Matrix, pPivot );
217 else
218 RealRowColElimination( Matrix, pPivot );
219 }
220 else
221 { ReorderingRequired = YES;
222 break; /* for loop */
223 }
224 }
225 if (NOT ReorderingRequired)
226 goto Done;
227 else
228 {
229 /*
230 * A pivot was not large enough to maintain accuracy,
231 * so a partial reordering is required.
232 */
233
234 #if ANNOTATE >= ON_STRANGE_BEHAVIOR
235 printf("Reordering, Step = %1d\n", Step);
236 #endif
237 }
238 } /* End of if(NOT Matrix->NeedsOrdering) */
239 else
240 {
241 /*
242 * This is the first time the matrix has been factored. These few statements
243 * indicate to the rest of the code that a full reodering is required rather
244 * than a partial reordering, which occurs during a failure of a fast
245 * factorization.
246 */
247 Step = 1;
248 if (NOT Matrix->RowsLinked)
249 spcLinkRows( Matrix );
250 if (NOT Matrix->InternalVectorsAllocated)
251 CreateInternalVectors( Matrix );
252 if (Matrix->Error >= spFATAL)
253 return Matrix->Error;
254 }
255
256 /* Form initial Markowitz products. */
257 CountMarkowitz( Matrix, RHS, Step );
258 MarkowitzProducts( Matrix, Step );
259 Matrix->MaxRowCountInLowerTri = -1;
260
261 /* Perform reordering and factorization. */
262 for (; Step <= Size; Step++)
263 { pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
264 if (pPivot == NULL) return MatrixIsSingular( Matrix, Step );
265 ExchangeRowsAndCols( Matrix, pPivot, Step );
266
267 if (Matrix->Complex)
268 ComplexRowColElimination( Matrix, pPivot );
269 else
270 RealRowColElimination( Matrix, pPivot );
271
272 if (Matrix->Error >= spFATAL) return Matrix->Error;
273 UpdateMarkowitzNumbers( Matrix, pPivot );
274
275 #if ANNOTATE == FULL
276 WriteStatus( Matrix, Step );
277 #endif
278 }
279
280 Done:
281 Matrix->NeedsOrdering = NO;
282 Matrix->Reordered = YES;
283 Matrix->Factored = YES;
284
285 return Matrix->Error;
286 }
287
288
289
290
291
292
293
294 /*
295 * FACTOR MATRIX
296 *
297 * This routine is the companion routine to spOrderAndFactor().
298 * Unlike spOrderAndFactor(), spFactor() cannot change the ordering.
299 * It is also faster than spOrderAndFactor(). The standard way of
300 * using these two routines is to first use spOrderAndFactor() for the
301 * initial factorization. For subsequent factorizations, spFactor()
302 * is used if there is some assurance that little growth will occur
303 * (say for example, that the matrix is diagonally dominant). If
304 * spFactor() is called for the initial factorization of the matrix,
305 * then spOrderAndFactor() is automatically called with the default
306 * threshold. This routine uses "row at a time" LU factorization.
307 * Pivots are associated with the lower triangular matrix and the
308 * diagonals of the upper triangular matrix are ones.
309 *
310 * >>> Returned:
311 * The error code is returned. Possible errors are listed below.
312 *
313 * >>> Arguments:
314 * Matrix <input> (char *)
315 * Pointer to matrix.
316 *
317 * >>> Possible errors:
318 * spNO_MEMORY
319 * spSINGULAR
320 * spZERO_DIAG
321 * spSMALL_PIVOT
322 * Error is cleared in this function.
323 */
324
325 int
spFactor(eMatrix)326 spFactor( eMatrix )
327
328 char *eMatrix;
329 {
330 MatrixPtr Matrix = (MatrixPtr)eMatrix;
331 register ElementPtr pElement;
332 register ElementPtr pColumn;
333 register int Step, Size;
334 RealNumber Mult;
335
336 /* Begin `spFactor'. */
337 ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
338
339 if (Matrix->NeedsOrdering)
340 { return spOrderAndFactor( eMatrix, (RealVector)NULL,
341 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT );
342 }
343 if (NOT Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION );
344 #if spCOMPLEX
345 if (Matrix->Complex) return FactorComplexMatrix( Matrix );
346 #endif
347
348 #if REAL
349 Size = Matrix->Size;
350
351 if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
352 Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;
353
354 /* Start factorization. */
355 for (Step = 2; Step <= Size; Step++)
356 { if (Matrix->DoRealDirect[Step])
357 { /* Update column using direct addressing scatter-gather. */
358 register RealNumber *Dest = (RealNumber *)Matrix->Intermediate;
359
360 /* Scatter. */
361 pElement = Matrix->FirstInCol[Step];
362 while (pElement != NULL)
363 { Dest[pElement->Row] = pElement->Real;
364 pElement = pElement->NextInCol;
365 }
366
367 /* Update column. */
368 pColumn = Matrix->FirstInCol[Step];
369 while (pColumn->Row < Step)
370 { pElement = Matrix->Diag[pColumn->Row];
371 pColumn->Real = Dest[pColumn->Row] * pElement->Real;
372 while ((pElement = pElement->NextInCol) != NULL)
373 Dest[pElement->Row] -= pColumn->Real * pElement->Real;
374 pColumn = pColumn->NextInCol;
375 }
376
377 /* Gather. */
378 pElement = Matrix->Diag[Step]->NextInCol;
379 while (pElement != NULL)
380 { pElement->Real = Dest[pElement->Row];
381 pElement = pElement->NextInCol;
382 }
383
384 /* Check for singular matrix. */
385 if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step );
386 Matrix->Diag[Step]->Real = 1.0 / Dest[Step];
387 }
388 else
389 { /* Update column using indirect addressing scatter-gather. */
390 register RealNumber **pDest = (RealNumber **)Matrix->Intermediate;
391
392 /* Scatter. */
393 pElement = Matrix->FirstInCol[Step];
394 while (pElement != NULL)
395 { pDest[pElement->Row] = &pElement->Real;
396 pElement = pElement->NextInCol;
397 }
398
399 /* Update column. */
400 pColumn = Matrix->FirstInCol[Step];
401 while (pColumn->Row < Step)
402 { pElement = Matrix->Diag[pColumn->Row];
403 Mult = (*pDest[pColumn->Row] *= pElement->Real);
404 while ((pElement = pElement->NextInCol) != NULL)
405 *pDest[pElement->Row] -= Mult * pElement->Real;
406 pColumn = pColumn->NextInCol;
407 }
408
409 /* Check for singular matrix. */
410 if (Matrix->Diag[Step]->Real == 0.0)
411 return ZeroPivot( Matrix, Step );
412 Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real;
413 }
414 }
415
416 Matrix->Factored = YES;
417 return (Matrix->Error = spOKAY);
418 #endif /* REAL */
419 }
420
421
422
423
424
425
426 #if spCOMPLEX
427 /*
428 * FACTOR COMPLEX MATRIX
429 *
430 * This routine is the companion routine to spFactor(), it
431 * handles complex matrices. It is otherwise identical.
432 *
433 * >>> Returned:
434 * The error code is returned. Possible errors are listed below.
435 *
436 * >>> Arguments:
437 * Matrix <input> (char *)
438 * Pointer to matrix.
439 *
440 * >>> Possible errors:
441 * spSINGULAR
442 * Error is cleared in this function.
443 */
444
445 static int
FactorComplexMatrix(Matrix)446 FactorComplexMatrix( Matrix )
447
448 MatrixPtr Matrix;
449 {
450 register ElementPtr pElement;
451 register ElementPtr pColumn;
452 register int Step, Size;
453 ComplexNumber Mult, Pivot;
454
455 /* Begin `FactorComplexMatrix'. */
456 ASSERT(Matrix->Complex);
457
458 Size = Matrix->Size;
459 pElement = Matrix->Diag[1];
460 if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 );
461 /* Cmplx expr: *pPivot = 1.0 / *pPivot. */
462 CMPLX_RECIPROCAL( *pElement, *pElement );
463
464 /* Start factorization. */
465 for (Step = 2; Step <= Size; Step++)
466 { if (Matrix->DoCmplxDirect[Step])
467 { /* Update column using direct addressing scatter-gather. */
468 register ComplexNumber *Dest;
469 Dest = (ComplexNumber *)Matrix->Intermediate;
470
471 /* Scatter. */
472 pElement = Matrix->FirstInCol[Step];
473 while (pElement != NULL)
474 { Dest[pElement->Row] = *(ComplexNumber *)pElement;
475 pElement = pElement->NextInCol;
476 }
477
478 /* Update column. */
479 pColumn = Matrix->FirstInCol[Step];
480 while (pColumn->Row < Step)
481 { pElement = Matrix->Diag[pColumn->Row];
482 /* Cmplx expr: Mult = Dest[pColumn->Row] * (1.0 / *pPivot). */
483 CMPLX_MULT(Mult, Dest[pColumn->Row], *pElement);
484 CMPLX_ASSIGN(*pColumn, Mult);
485 while ((pElement = pElement->NextInCol) != NULL)
486 { /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */
487 CMPLX_MULT_SUBT_ASSIGN(Dest[pElement->Row],Mult,*pElement);
488 }
489 pColumn = pColumn->NextInCol;
490 }
491
492 /* Gather. */
493 pElement = Matrix->Diag[Step]->NextInCol;
494 while (pElement != NULL)
495 { *(ComplexNumber *)pElement = Dest[pElement->Row];
496 pElement = pElement->NextInCol;
497 }
498
499 /* Check for singular matrix. */
500 Pivot = Dest[Step];
501 if (CMPLX_1_NORM(Pivot) == 0.0) return ZeroPivot( Matrix, Step );
502 CMPLX_RECIPROCAL( *Matrix->Diag[Step], Pivot );
503 }
504 else
505 { /* Update column using direct addressing scatter-gather. */
506 register ComplexNumber **pDest;
507 pDest = (ComplexNumber **)Matrix->Intermediate;
508
509 /* Scatter. */
510 pElement = Matrix->FirstInCol[Step];
511 while (pElement != NULL)
512 { pDest[pElement->Row] = (ComplexNumber *)pElement;
513 pElement = pElement->NextInCol;
514 }
515
516 /* Update column. */
517 pColumn = Matrix->FirstInCol[Step];
518 while (pColumn->Row < Step)
519 { pElement = Matrix->Diag[pColumn->Row];
520 /* Cmplx expr: Mult = *pDest[pColumn->Row] * (1.0 / *pPivot). */
521 CMPLX_MULT(Mult, *pDest[pColumn->Row], *pElement);
522 CMPLX_ASSIGN(*pDest[pColumn->Row], Mult);
523 while ((pElement = pElement->NextInCol) != NULL)
524 { /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */
525 CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row],Mult,*pElement);
526 }
527 pColumn = pColumn->NextInCol;
528 }
529
530 /* Check for singular matrix. */
531 pElement = Matrix->Diag[Step];
532 if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, Step );
533 CMPLX_RECIPROCAL( *pElement, *pElement );
534 }
535 }
536
537 Matrix->Factored = YES;
538 return (Matrix->Error = spOKAY);
539 }
540 #endif /* spCOMPLEX */
541
542
543
544
545
546
547 /*
548 * PARTITION MATRIX
549 *
550 * This routine determines the cost to factor each row using both
551 * direct and indirect addressing and decides, on a row-by-row basis,
552 * which addressing mode is fastest. This information is used in
553 * spFactor() to speed the factorization.
554 *
555 * When factoring a previously ordered matrix using spFactor(), Sparse
556 * operates on a row-at-a-time basis. For speed, on each step, the
557 * row being updated is copied into a full vector and the operations
558 * are performed on that vector. This can be done one of two ways,
559 * either using direct addressing or indirect addressing. Direct
560 * addressing is fastest when the matrix is relatively dense and
561 * indirect addressing is best when the matrix is quite sparse. The
562 * user selects the type of partition used with Mode. If Mode is set
563 * to spDIRECT_PARTITION, then the all rows are placed in the direct
564 * addressing partition. Similarly, if Mode is set to
565 * spINDIRECT_PARTITION, then the all rows are placed in the indirect
566 * addressing partition. By setting Mode to spAUTO_PARTITION, the
567 * user allows Sparse to select the partition for each row
568 * individually. spFactor() generally runs faster if Sparse is
569 * allowed to choose its own partitioning, however choosing a
570 * partition is expensive. The time required to choose a partition is
571 * of the same order of the cost to factor the matrix. If you plan to
572 * factor a large number of matrices with the same structure, it is
573 * best to let Sparse choose the partition. Otherwise, you should
574 * choose the partition based on the predicted density of the matrix.
575 *
576 * >>> Arguments:
577 * Matrix <input> (char *)
578 * Pointer to matrix.
579 * Mode <input> (int)
580 * Mode must be one of three special codes: spDIRECT_PARTITION,
581 * spINDIRECT_PARTITION, or spAUTO_PARTITION.
582 */
583
584 void
spPartition(eMatrix,Mode)585 spPartition( eMatrix, Mode )
586
587 char *eMatrix;
588 int Mode;
589 {
590 MatrixPtr Matrix = (MatrixPtr)eMatrix;
591 register ElementPtr pElement, pColumn;
592 register int Step, Size;
593 register int *Nc, *No, *Nm;
594 BOOLEAN *DoRealDirect, *DoCmplxDirect;
595
596 /* Begin `spPartition'. */
597 ASSERT( IS_SPARSE( Matrix ) );
598 if (Matrix->Partitioned) return;
599 Size = Matrix->Size;
600 DoRealDirect = Matrix->DoRealDirect;
601 DoCmplxDirect = Matrix->DoCmplxDirect;
602 Matrix->Partitioned = YES;
603
604 /* If partition is specified by the user, this is easy. */
605 if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION;
606 if (Mode == spDIRECT_PARTITION)
607 { for (Step = 1; Step <= Size; Step++)
608 #if REAL
609 DoRealDirect[Step] = YES;
610 #endif
611 #if spCOMPLEX
612 DoCmplxDirect[Step] = YES;
613 #endif
614 return;
615 }
616 else if (Mode == spINDIRECT_PARTITION)
617 { for (Step = 1; Step <= Size; Step++)
618 #if REAL
619 DoRealDirect[Step] = NO;
620 #endif
621 #if spCOMPLEX
622 DoCmplxDirect[Step] = NO;
623 #endif
624 return;
625 }
626 else ASSERT( Mode == spAUTO_PARTITION );
627
628 /* Otherwise, count all operations needed in when factoring matrix. */
629 Nc = (int *)Matrix->MarkowitzRow;
630 No = (int *)Matrix->MarkowitzCol;
631 Nm = (int *)Matrix->MarkowitzProd;
632
633 /* Start mock-factorization. */
634 for (Step = 1; Step <= Size; Step++)
635 { Nc[Step] = No[Step] = Nm[Step] = 0;
636
637 pElement = Matrix->FirstInCol[Step];
638 while (pElement != NULL)
639 { Nc[Step]++;
640 pElement = pElement->NextInCol;
641 }
642
643 pColumn = Matrix->FirstInCol[Step];
644 while (pColumn->Row < Step)
645 { pElement = Matrix->Diag[pColumn->Row];
646 Nm[Step]++;
647 while ((pElement = pElement->NextInCol) != NULL)
648 No[Step]++;
649 pColumn = pColumn->NextInCol;
650 }
651 }
652
653 for (Step = 1; Step <= Size; Step++)
654 {
655 /*
656 * The following are just estimates based on a count on the number of
657 * machine instructions used on each machine to perform the various
658 * tasks. It was assumed that each machine instruction required the
659 * same amount of time (I don't believe this is true for the VAX, and
660 * have no idea if this is true for the 68000 family). For optimum
661 * performance, these numbers should be tuned to the machine.
662 * Nc is the number of nonzero elements in the column.
663 * Nm is the number of multipliers in the column.
664 * No is the number of operations in the inner loop.
665 */
666
667 #define generic
668 #ifdef hp9000s300
669 #if REAL
670 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
671 #endif
672 #if spCOMPLEX
673 /* On the hp350, it is never profitable to use direct for complex. */
674 DoCmplxDirect[Step] = NO;
675 #endif
676 #undef generic
677 #endif
678
679 #ifdef vax
680 #if REAL
681 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
682 #endif
683 #if spCOMPLEX
684 DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
685 #endif
686 #undef generic
687 #endif
688
689 #ifdef generic
690 #if REAL
691 DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
692 #endif
693 #if spCOMPLEX
694 DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
695 #endif
696 #undef generic
697 #endif
698 }
699
700 #if (ANNOTATE == FULL)
701 { int Ops = 0;
702 for (Step = 1; Step <= Size; Step++)
703 Ops += No[Step];
704 printf("Operation count for inner loop of factorization = %d.\n", Ops);
705 }
706 #endif
707 return;
708 }
709
710
711
712
713
714
715
716 /*
717 * CREATE INTERNAL VECTORS
718 *
719 * Creates the Markowitz and Intermediate vectors.
720 *
721 * >>> Arguments:
722 * Matrix <input> (MatrixPtr)
723 * Pointer to matrix.
724 *
725 * >>> Local variables:
726 * SizePlusOne (unsigned)
727 * Size of the arrays to be allocated.
728 *
729 * >>> Possible errors:
730 * spNO_MEMORY
731 */
732
733 static
CreateInternalVectors(Matrix)734 CreateInternalVectors( Matrix )
735
736 MatrixPtr Matrix;
737 {
738 int Size;
739
740 /* Begin `CreateInternalVectors'. */
741 /* Create Markowitz arrays. */
742 Size= Matrix->Size;
743
744 if (Matrix->MarkowitzRow == NULL)
745 { if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL)
746 Matrix->Error = spNO_MEMORY;
747 }
748 if (Matrix->MarkowitzCol == NULL)
749 { if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL)
750 Matrix->Error = spNO_MEMORY;
751 }
752 if (Matrix->MarkowitzProd == NULL)
753 { if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL)
754 Matrix->Error = spNO_MEMORY;
755 }
756
757 /* Create DoDirect vectors for use in spFactor(). */
758 #if REAL
759 if (Matrix->DoRealDirect == NULL)
760 { if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
761 Matrix->Error = spNO_MEMORY;
762 }
763 #endif
764 #if spCOMPLEX
765 if (Matrix->DoCmplxDirect == NULL)
766 { if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
767 Matrix->Error = spNO_MEMORY;
768 }
769 #endif
770
771 /* Create Intermediate vectors for use in MatrixSolve. */
772 #if spCOMPLEX
773 if (Matrix->Intermediate == NULL)
774 { if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL)
775 Matrix->Error = spNO_MEMORY;
776 }
777 #else
778 if (Matrix->Intermediate == NULL)
779 { if ((Matrix->Intermediate = ALLOC(RealNumber, Size+1)) == NULL)
780 Matrix->Error = spNO_MEMORY;
781 }
782 #endif
783
784 if (Matrix->Error != spNO_MEMORY)
785 Matrix->InternalVectorsAllocated = YES;
786 return 0;
787 }
788
789
790
791
792
793
794
795 /*
796 * COUNT MARKOWITZ
797 *
798 * Scans Matrix to determine the Markowitz counts for each row and column.
799 *
800 * >>> Arguments:
801 * Matrix <input> (MatrixPtr)
802 * Pointer to matrix.
803 * RHS <input> (RealVector)
804 * Representative right-hand side vector that is used to determine
805 * pivoting order when the right hand side vector is sparse. If
806 * RHS is a NULL pointer then the RHS vector is assumed to be full
807 * and it is not used when determining the pivoting order.
808 * Step <input> (int)
809 * Index of the diagonal currently being eliminated.
810 *
811 * >>> Local variables:
812 * Count (int)
813 * Temporary counting variable.
814 * ExtRow (int)
815 * The external row number that corresponds to I.
816 * pElement (ElementPtr)
817 * Pointer to matrix elements.
818 * Size (int)
819 * The size of the matrix.
820 */
821
822 static
CountMarkowitz(Matrix,RHS,Step)823 CountMarkowitz( Matrix, RHS, Step )
824
825 MatrixPtr Matrix;
826 register RealVector RHS;
827 int Step;
828 {
829 register int Count, I, Size = Matrix->Size;
830 register ElementPtr pElement;
831 int ExtRow;
832
833 /* Begin `CountMarkowitz'. */
834
835 /* Correct array pointer for ARRAY_OFFSET. */
836 #if NOT ARRAY_OFFSET
837 #if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX
838 if (RHS != NULL) --RHS;
839 #else
840 if (RHS != NULL)
841 { if (Matrix->Complex) RHS -= 2;
842 else --RHS;
843 }
844 #endif
845 #endif
846
847 /* Generate MarkowitzRow Count for each row. */
848 for (I = Step; I <= Size; I++)
849 {
850 /* Set Count to -1 initially to remove count due to pivot element. */
851 Count = -1;
852 pElement = Matrix->FirstInRow[I];
853 while (pElement != NULL AND pElement->Col < Step)
854 pElement = pElement->NextInRow;
855 while (pElement != NULL)
856 { Count++;
857 pElement = pElement->NextInRow;
858 }
859
860 /* Include nonzero elements in the RHS vector. */
861 ExtRow = Matrix->IntToExtRowMap[I];
862
863 #if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX
864 if (RHS != NULL)
865 if (RHS[ExtRow] != 0.0) Count++;
866 #else
867 if (RHS != NULL)
868 { if (Matrix->Complex)
869 { if ((RHS[2*ExtRow] != 0.0) OR (RHS[2*ExtRow+1] != 0.0))
870 Count++;
871 }
872 else if (RHS[I] != 0.0) Count++;
873 }
874 #endif
875 Matrix->MarkowitzRow[I] = Count;
876 }
877
878 /* Generate the MarkowitzCol count for each column. */
879 for (I = Step; I <= Size; I++)
880 {
881 /* Set Count to -1 initially to remove count due to pivot element. */
882 Count = -1;
883 pElement = Matrix->FirstInCol[I];
884 while (pElement != NULL AND pElement->Row < Step)
885 pElement = pElement->NextInCol;
886 while (pElement != NULL)
887 { Count++;
888 pElement = pElement->NextInCol;
889 }
890 Matrix->MarkowitzCol[I] = Count;
891 }
892 return 0;
893 }
894
895
896
897
898
899
900
901
902
903
904 /*
905 * MARKOWITZ PRODUCTS
906 *
907 * Calculates MarkowitzProduct for each diagonal element from the Markowitz
908 * counts.
909 *
910 * >>> Arguments:
911 * Matrix <input> (MatrixPtr)
912 * Pointer to matrix.
913 * Step <input> (int)
914 * Index of the diagonal currently being eliminated.
915 *
916 * >>> Local Variables:
917 * pMarkowitzProduct (long *)
918 * Pointer that points into MarkowitzProduct array. Is used to
919 * sequentially access entries quickly.
920 * pMarkowitzRow (int *)
921 * Pointer that points into MarkowitzRow array. Is used to sequentially
922 * access entries quickly.
923 * pMarkowitzCol (int *)
924 * Pointer that points into MarkowitzCol array. Is used to sequentially
925 * access entries quickly.
926 * Product (long)
927 * Temporary storage for Markowitz product./
928 * Size (int)
929 * The size of the matrix.
930 */
931
932 static
MarkowitzProducts(Matrix,Step)933 MarkowitzProducts( Matrix, Step )
934
935 MatrixPtr Matrix;
936 int Step;
937 {
938 register int I, *pMarkowitzRow, *pMarkowitzCol;
939 register long Product, *pMarkowitzProduct;
940 register int Size = Matrix->Size;
941 double fProduct;
942
943 /* Begin `MarkowitzProducts'. */
944 Matrix->Singletons = 0;
945
946 pMarkowitzProduct = &(Matrix->MarkowitzProd[Step]);
947 pMarkowitzRow = &(Matrix->MarkowitzRow[Step]);
948 pMarkowitzCol = &(Matrix->MarkowitzCol[Step]);
949
950 for (I = Step; I <= Size; I++)
951 {
952 /* If chance of overflow, use real numbers. */
953 if ((*pMarkowitzRow > LARGEST_SHORT_INTEGER AND *pMarkowitzCol != 0) OR
954 (*pMarkowitzCol > LARGEST_SHORT_INTEGER AND *pMarkowitzRow != 0))
955 { fProduct = (double)(*pMarkowitzRow++) * (double)(*pMarkowitzCol++);
956 if (fProduct >= LARGEST_LONG_INTEGER)
957 *pMarkowitzProduct++ = LARGEST_LONG_INTEGER;
958 else
959 *pMarkowitzProduct++ = fProduct;
960 }
961 else
962 { Product = *pMarkowitzRow++ * *pMarkowitzCol++;
963 if ((*pMarkowitzProduct++ = Product) == 0)
964 Matrix->Singletons++;
965 }
966 }
967 return 0;
968 }
969
970
971
972
973
974
975
976
977
978
979
980 /*
981 * SEARCH FOR BEST PIVOT
982 *
983 * Performs a search to determine the element with the lowest Markowitz
984 * Product that is also acceptable. An acceptable element is one that is
985 * larger than the AbsThreshold and at least as large as RelThreshold times
986 * the largest element in the same column. The first step is to look for
987 * singletons if any exist. If none are found, then all the diagonals are
988 * searched. The diagonal is searched once quickly using the assumption that
989 * elements on the diagonal are large compared to other elements in their
990 * column, and so the pivot can be chosen only on the basis of the Markowitz
991 * criterion. After a element has been chosen to be pivot on the basis of
992 * its Markowitz product, it is checked to see if it is large enough.
993 * Waiting to the end of the Markowitz search to check the size of a pivot
994 * candidate saves considerable time, but is not guaranteed to find an
995 * acceptable pivot. Thus if unsuccessful a second pass of the diagonal is
996 * made. This second pass checks to see if an element is large enough during
997 * the search, not after it. If still no acceptable pivot candidate has
998 * been found, the search expands to cover the entire matrix.
999 *
1000 * >>> Returned:
1001 * A pointer to the element chosen to be pivot. If every element in the
1002 * matrix is zero, then NULL is returned.
1003 *
1004 * >>> Arguments:
1005 * Matrix <input> (MatrixPtr)
1006 * Pointer to matrix.
1007 * Step <input> (int)
1008 * The row and column number of the beginning of the reduced submatrix.
1009 *
1010 * >>> Local variables:
1011 * ChosenPivot (ElementPtr)
1012 * Pointer to element that has been chosen to be the pivot.
1013 *
1014 * >>> Possible errors:
1015 * spSINGULAR
1016 * spSMALL_PIVOT
1017 */
1018
1019 static ElementPtr
SearchForPivot(Matrix,Step,DiagPivoting)1020 SearchForPivot( Matrix, Step, DiagPivoting )
1021
1022 MatrixPtr Matrix;
1023 int Step, DiagPivoting;
1024 {
1025 register ElementPtr ChosenPivot;
1026 ElementPtr SearchForSingleton();
1027 ElementPtr QuicklySearchDiagonal();
1028 ElementPtr SearchDiagonal();
1029 ElementPtr SearchEntireMatrix();
1030
1031 /* Begin `SearchForPivot'. */
1032
1033 /* If singletons exist, look for an acceptable one to use as pivot. */
1034 if (Matrix->Singletons)
1035 { ChosenPivot = SearchForSingleton( Matrix, Step );
1036 if (ChosenPivot != NULL)
1037 { Matrix->PivotSelectionMethod = 's';
1038 return ChosenPivot;
1039 }
1040 }
1041
1042 #if DIAGONAL_PIVOTING
1043 if (DiagPivoting)
1044 {
1045 /*
1046 * Either no singletons exist or they weren't acceptable. Take quick first
1047 * pass at searching diagonal. First search for element on diagonal of
1048 * remaining submatrix with smallest Markowitz product, then check to see
1049 * if it okay numerically. If not, QuicklySearchDiagonal fails.
1050 */
1051 ChosenPivot = QuicklySearchDiagonal( Matrix, Step );
1052 if (ChosenPivot != NULL)
1053 { Matrix->PivotSelectionMethod = 'q';
1054 return ChosenPivot;
1055 }
1056
1057 /*
1058 * Quick search of diagonal failed, carefully search diagonal and check each
1059 * pivot candidate numerically before even tentatively accepting it.
1060 */
1061 ChosenPivot = SearchDiagonal( Matrix, Step );
1062 if (ChosenPivot != NULL)
1063 { Matrix->PivotSelectionMethod = 'd';
1064 return ChosenPivot;
1065 }
1066 }
1067 #endif /* DIAGONAL_PIVOTING */
1068
1069 /* No acceptable pivot found yet, search entire matrix. */
1070 ChosenPivot = SearchEntireMatrix( Matrix, Step );
1071 Matrix->PivotSelectionMethod = 'e';
1072
1073 return ChosenPivot;
1074 }
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084 /*
1085 * SEARCH FOR SINGLETON TO USE AS PIVOT
1086 *
1087 * Performs a search to find a singleton to use as the pivot. The
1088 * first acceptable singleton is used. A singleton is acceptable if
1089 * it is larger in magnitude than the AbsThreshold and larger
1090 * than RelThreshold times the largest of any other elements in the same
1091 * column. It may seem that a singleton need not satisfy the
1092 * relative threshold criterion, however it is necessary to prevent
1093 * excessive growth in the RHS from resulting in overflow during the
1094 * forward and backward substitution. A singleton does not need to
1095 * be on the diagonal to be selected.
1096 *
1097 * >>> Returned:
1098 * A pointer to the singleton chosen to be pivot. In no singleton is
1099 * acceptable, return NULL.
1100 *
1101 * >>> Arguments:
1102 * Matrix <input> (MatrixPtr)
1103 * Pointer to matrix.
1104 * Step <input> (int)
1105 * Index of the diagonal currently being eliminated.
1106 *
1107 * >>> Local variables:
1108 * ChosenPivot (ElementPtr)
1109 * Pointer to element that has been chosen to be the pivot.
1110 * PivotMag (RealNumber)
1111 * Magnitude of ChosenPivot.
1112 * Singletons (int)
1113 * The count of the number of singletons that can be used as pivots.
1114 * A local version of Matrix->Singletons.
1115 * pMarkowitzProduct (long *)
1116 * Pointer that points into MarkowitzProduct array. It is used to quickly
1117 * access successive Markowitz products.
1118 */
1119
1120 static ElementPtr
SearchForSingleton(Matrix,Step)1121 SearchForSingleton( Matrix, Step )
1122
1123 MatrixPtr Matrix;
1124 int Step;
1125 {
1126 register ElementPtr ChosenPivot;
1127 register int I;
1128 register long *pMarkowitzProduct;
1129 int Singletons;
1130 RealNumber PivotMag, FindBiggestInColExclude();
1131
1132 /* Begin `SearchForSingleton'. */
1133 /* Initialize pointer that is to scan through MarkowitzProduct vector. */
1134 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+1]);
1135 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
1136
1137 /* Decrement the count of available singletons, on the assumption that an
1138 * acceptable one will be found. */
1139 Singletons = Matrix->Singletons--;
1140
1141 /*
1142 * Assure that following while loop will always terminate, this is just
1143 * preventive medicine, if things are working right this should never
1144 * be needed.
1145 */
1146 Matrix->MarkowitzProd[Step-1] = 0;
1147
1148 while (Singletons-- > 0)
1149 {
1150 /* Singletons exist, find them. */
1151
1152 /*
1153 * This is tricky. Am using a pointer to sequentially step through the
1154 * MarkowitzProduct array. Search terminates when singleton (Product = 0)
1155 * is found. Note that the conditional in the while statement
1156 * ( *pMarkowitzProduct ) is true as long as the MarkowitzProduct is not
1157 * equal to zero. The row (and column) index on the diagonal is then
1158 * calculated by subtracting the pointer to the Markowitz product of
1159 * the first diagonal from the pointer to the Markowitz product of the
1160 * desired element, the singleton.
1161 *
1162 * Search proceeds from the end (high row and column numbers) to the
1163 * beginning (low row and column numbers) so that rows and columns with
1164 * large Markowitz products will tend to be move to the bottom of the
1165 * matrix. However, choosing Diag[Step] is desirable because it would
1166 * require no row and column interchanges, so inspect it first by
1167 * putting its Markowitz product at the end of the MarkowitzProd
1168 * vector.
1169 */
1170
1171 while ( *pMarkowitzProduct-- )
1172 { /*
1173 * N bottles of beer on the wall;
1174 * N bottles of beer.
1175 * you take one down and pass it around;
1176 * N-1 bottles of beer on the wall.
1177 */
1178 }
1179 I = pMarkowitzProduct - Matrix->MarkowitzProd + 1;
1180
1181 /* Assure that I is valid. */
1182 if (I < Step) break; /* while (Singletons-- > 0) */
1183 if (I > Matrix->Size) I = Step;
1184
1185 /* Singleton has been found in either/both row or/and column I. */
1186 if ((ChosenPivot = Matrix->Diag[I]) != NULL)
1187 {
1188 /* Singleton lies on the diagonal. */
1189 PivotMag = ELEMENT_MAG(ChosenPivot);
1190 if
1191 ( PivotMag > Matrix->AbsThreshold AND
1192 PivotMag > Matrix->RelThreshold *
1193 FindBiggestInColExclude( Matrix, ChosenPivot, Step )
1194 ) return ChosenPivot;
1195 }
1196 else
1197 {
1198 /* Singleton does not lie on diagonal, find it. */
1199 if (Matrix->MarkowitzCol[I] == 0)
1200 { ChosenPivot = Matrix->FirstInCol[I];
1201 while ((ChosenPivot != NULL) AND (ChosenPivot->Row < Step))
1202 ChosenPivot = ChosenPivot->NextInCol;
1203 PivotMag = ELEMENT_MAG(ChosenPivot);
1204 if
1205 ( PivotMag > Matrix->AbsThreshold AND
1206 PivotMag > Matrix->RelThreshold *
1207 FindBiggestInColExclude( Matrix, ChosenPivot,
1208 Step )
1209 ) return ChosenPivot;
1210 else
1211 { if (Matrix->MarkowitzRow[I] == 0)
1212 { ChosenPivot = Matrix->FirstInRow[I];
1213 while((ChosenPivot != NULL) AND (ChosenPivot->Col<Step))
1214 ChosenPivot = ChosenPivot->NextInRow;
1215 PivotMag = ELEMENT_MAG(ChosenPivot);
1216 if
1217 ( PivotMag > Matrix->AbsThreshold AND
1218 PivotMag > Matrix->RelThreshold *
1219 FindBiggestInColExclude( Matrix,
1220 ChosenPivot,
1221 Step )
1222 ) return ChosenPivot;
1223 }
1224 }
1225 }
1226 else
1227 { ChosenPivot = Matrix->FirstInRow[I];
1228 while ((ChosenPivot != NULL) AND (ChosenPivot->Col < Step))
1229 ChosenPivot = ChosenPivot->NextInRow;
1230 PivotMag = ELEMENT_MAG(ChosenPivot);
1231 if
1232 ( PivotMag > Matrix->AbsThreshold AND
1233 PivotMag > Matrix->RelThreshold *
1234 FindBiggestInColExclude( Matrix, ChosenPivot,
1235 Step )
1236 ) return ChosenPivot;
1237 }
1238 }
1239 /* Singleton not acceptable (too small), try another. */
1240 } /* end of while(lSingletons>0) */
1241
1242 /*
1243 * All singletons were unacceptable. Restore Matrix->Singletons count.
1244 * Initial assumption that an acceptable singleton would be found was wrong.
1245 */
1246 Matrix->Singletons++;
1247 return NULL;
1248 }
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261 #if DIAGONAL_PIVOTING
1262 #if MODIFIED_MARKOWITZ
1263 /*
1264 * QUICK SEARCH OF DIAGONAL FOR PIVOT WITH MODIFIED MARKOWITZ CRITERION
1265 *
1266 * Searches the diagonal looking for the best pivot. For a pivot to be
1267 * acceptable it must be larger than the pivot RelThreshold times the largest
1268 * element in its reduced column. Among the acceptable diagonals, the
1269 * one with the smallest MarkowitzProduct is sought. Search terminates
1270 * early if a diagonal is found with a MarkowitzProduct of one and its
1271 * magnitude is larger than the other elements in its row and column.
1272 * Since its MarkowitzProduct is one, there is only one other element in
1273 * both its row and column, and, as a condition for early termination,
1274 * these elements must be located symmetricly in the matrix. If a tie
1275 * occurs between elements of equal MarkowitzProduct, then the element with
1276 * the largest ratio between its magnitude and the largest element in its
1277 * column is used. The search will be terminated after a given number of
1278 * ties have occurred and the best (largest ratio) of the tied element will
1279 * be used as the pivot. The number of ties that will trigger an early
1280 * termination is MinMarkowitzProduct * TIES_MULTIPLIER.
1281 *
1282 * >>> Returned:
1283 * A pointer to the diagonal element chosen to be pivot. If no diagonal is
1284 * acceptable, a NULL is returned.
1285 *
1286 * >>> Arguments:
1287 * Step <input> (int)
1288 * Index of the diagonal currently being eliminated.
1289 *
1290 * >>> Local variables:
1291 * ChosenPivot (ElementPtr)
1292 * Pointer to the element that has been chosen to be the pivot.
1293 * LargestOffDiagonal (RealNumber)
1294 * Magnitude of the largest of the off-diagonal terms associated with
1295 * a diagonal with MarkowitzProduct equal to one.
1296 * Magnitude (RealNumber)
1297 * Absolute value of diagonal element.
1298 * MaxRatio (RealNumber)
1299 * Among the elements tied with the smallest Markowitz product, MaxRatio
1300 * is the best (smallest) ratio of LargestInCol to the diagonal Magnitude
1301 * found so far. The smaller the ratio, the better numerically the
1302 * element will be as pivot.
1303 * MinMarkowitzProduct (long)
1304 * Smallest Markowitz product found of pivot candidates that lie along
1305 * diagonal.
1306 * NumberOfTies (int)
1307 * A count of the number of Markowitz ties that have occurred at current
1308 * MarkowitzProduct.
1309 * pDiag (ElementPtr)
1310 * Pointer to current diagonal element.
1311 * pMarkowitzProduct (long *)
1312 * Pointer that points into MarkowitzProduct array. It is used to quickly
1313 * access successive Markowitz products.
1314 * Ratio (RealNumber)
1315 * For the current pivot candidate, Ratio is the ratio of the largest
1316 * element in its column (excluding itself) to its magnitude.
1317 * TiedElements (ElementPtr[])
1318 * Array of pointers to the elements with the minimum Markowitz
1319 * product.
1320 * pOtherInCol (ElementPtr)
1321 * When there is only one other element in a column other than the
1322 * diagonal, pOtherInCol is used to point to it. Used when Markowitz
1323 * product is to determine if off diagonals are placed symmetricly.
1324 * pOtherInRow (ElementPtr)
1325 * When there is only one other element in a row other than the diagonal,
1326 * pOtherInRow is used to point to it. Used when Markowitz product is
1327 * to determine if off diagonals are placed symmetricly.
1328 */
1329
1330 static ElementPtr
QuicklySearchDiagonal(Matrix,Step)1331 QuicklySearchDiagonal( Matrix, Step )
1332
1333 MatrixPtr Matrix;
1334 int Step;
1335 {
1336 register long MinMarkowitzProduct, *pMarkowitzProduct;
1337 register ElementPtr pDiag, pOtherInRow, pOtherInCol;
1338 int I, NumberOfTies;
1339 ElementPtr ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1];
1340 RealNumber Magnitude, LargestInCol, Ratio, MaxRatio;
1341 RealNumber LargestOffDiagonal;
1342 RealNumber FindBiggestInColExclude();
1343
1344 /* Begin `QuicklySearchDiagonal'. */
1345 NumberOfTies = -1;
1346 MinMarkowitzProduct = LARGEST_LONG_INTEGER;
1347 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
1348 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
1349
1350 /* Assure that following while loop will always terminate. */
1351 Matrix->MarkowitzProd[Step-1] = -1;
1352
1353 /*
1354 * This is tricky. Am using a pointer in the inner while loop to
1355 * sequentially step through the MarkowitzProduct array. Search
1356 * terminates when the Markowitz product of zero placed at location
1357 * Step-1 is found. The row (and column) index on the diagonal is then
1358 * calculated by subtracting the pointer to the Markowitz product of
1359 * the first diagonal from the pointer to the Markowitz product of the
1360 * desired element. The outer for loop is infinite, broken by using
1361 * break.
1362 *
1363 * Search proceeds from the end (high row and column numbers) to the
1364 * beginning (low row and column numbers) so that rows and columns with
1365 * large Markowitz products will tend to be move to the bottom of the
1366 * matrix. However, choosing Diag[Step] is desirable because it would
1367 * require no row and column interchanges, so inspect it first by
1368 * putting its Markowitz product at the end of the MarkowitzProd
1369 * vector.
1370 */
1371
1372 for(;;) /* Endless for loop. */
1373 { while (MinMarkowitzProduct < *(--pMarkowitzProduct))
1374 { /*
1375 * N bottles of beer on the wall;
1376 * N bottles of beer.
1377 * You take one down and pass it around;
1378 * N-1 bottles of beer on the wall.
1379 */
1380 }
1381
1382 I = pMarkowitzProduct - Matrix->MarkowitzProd;
1383
1384 /* Assure that I is valid; if I < Step, terminate search. */
1385 if (I < Step) break; /* Endless for loop */
1386 if (I > Matrix->Size) I = Step;
1387
1388 if ((pDiag = Matrix->Diag[I]) == NULL)
1389 continue; /* Endless for loop */
1390 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
1391 continue; /* Endless for loop */
1392
1393 if (*pMarkowitzProduct == 1)
1394 {
1395 /* Case where only one element exists in row and column other than diagonal. */
1396
1397 /* Find off diagonal elements. */
1398 pOtherInRow = pDiag->NextInRow;
1399 pOtherInCol = pDiag->NextInCol;
1400 if (pOtherInRow == NULL AND pOtherInCol == NULL)
1401 { pOtherInRow = Matrix->FirstInRow[I];
1402 while(pOtherInRow != NULL)
1403 { if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I)
1404 break;
1405 pOtherInRow = pOtherInRow->NextInRow;
1406 }
1407 pOtherInCol = Matrix->FirstInCol[I];
1408 while(pOtherInCol != NULL)
1409 { if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I)
1410 break;
1411 pOtherInCol = pOtherInCol->NextInCol;
1412 }
1413 }
1414
1415 /* Accept diagonal as pivot if diagonal is larger than off diagonals and the
1416 * off diagonals are placed symmetricly. */
1417 if (pOtherInRow != NULL AND pOtherInCol != NULL)
1418 { if (pOtherInRow->Col == pOtherInCol->Row)
1419 { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
1420 ELEMENT_MAG(pOtherInCol));
1421 if (Magnitude >= LargestOffDiagonal)
1422 {
1423 /* Accept pivot, it is unlikely to contribute excess error. */
1424 return pDiag;
1425 }
1426 }
1427 }
1428 }
1429
1430 if (*pMarkowitzProduct < MinMarkowitzProduct)
1431 {
1432 /* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
1433 TiedElements[0] = pDiag;
1434 MinMarkowitzProduct = *pMarkowitzProduct;
1435 NumberOfTies = 0;
1436 }
1437 else
1438 {
1439 /* This case handles Markowitz ties. */
1440 if (NumberOfTies < MAX_MARKOWITZ_TIES)
1441 { TiedElements[++NumberOfTies] = pDiag;
1442 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
1443 break; /* Endless for loop */
1444 }
1445 }
1446 } /* End of endless for loop. */
1447
1448 /* Test to see if any element was chosen as a pivot candidate. */
1449 if (NumberOfTies < 0)
1450 return NULL;
1451
1452 /* Determine which of tied elements is best numerically. */
1453 ChosenPivot = NULL;
1454 MaxRatio = 1.0 / Matrix->RelThreshold;
1455
1456 for (I = 0; I <= NumberOfTies; I++)
1457 { pDiag = TiedElements[I];
1458 Magnitude = ELEMENT_MAG(pDiag);
1459 LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
1460 Ratio = LargestInCol / Magnitude;
1461 if (Ratio < MaxRatio)
1462 { ChosenPivot = pDiag;
1463 MaxRatio = Ratio;
1464 }
1465 }
1466 return ChosenPivot;
1467 }
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478 #else /* Not MODIFIED_MARKOWITZ */
1479 /*
1480 * QUICK SEARCH OF DIAGONAL FOR PIVOT WITH CONVENTIONAL MARKOWITZ
1481 * CRITERION
1482 *
1483 * Searches the diagonal looking for the best pivot. For a pivot to be
1484 * acceptable it must be larger than the pivot RelThreshold times the largest
1485 * element in its reduced column. Among the acceptable diagonals, the
1486 * one with the smallest MarkowitzProduct is sought. Search terminates
1487 * early if a diagonal is found with a MarkowitzProduct of one and its
1488 * magnitude is larger than the other elements in its row and column.
1489 * Since its MarkowitzProduct is one, there is only one other element in
1490 * both its row and column, and, as a condition for early termination,
1491 * these elements must be located symmetricly in the matrix.
1492 *
1493 * >>> Returned:
1494 * A pointer to the diagonal element chosen to be pivot. If no diagonal is
1495 * acceptable, a NULL is returned.
1496 *
1497 * >>> Arguments:
1498 * Matrix <input> (MatrixPtr)
1499 * Pointer to matrix.
1500 * Step <input> (int)
1501 * Index of the diagonal currently being eliminated.
1502 *
1503 * >>> Local variables:
1504 * ChosenPivot (ElementPtr)
1505 * Pointer to the element that has been chosen to be the pivot.
1506 * LargestOffDiagonal (RealNumber)
1507 * Magnitude of the largest of the off-diagonal terms associated with
1508 * a diagonal with MarkowitzProduct equal to one.
1509 * Magnitude (RealNumber)
1510 * Absolute value of diagonal element.
1511 * MinMarkowitzProduct (long)
1512 * Smallest Markowitz product found of pivot candidates which are
1513 * acceptable.
1514 * pDiag (ElementPtr)
1515 * Pointer to current diagonal element.
1516 * pMarkowitzProduct (long *)
1517 * Pointer that points into MarkowitzProduct array. It is used to quickly
1518 * access successive Markowitz products.
1519 * pOtherInCol (ElementPtr)
1520 * When there is only one other element in a column other than the
1521 * diagonal, pOtherInCol is used to point to it. Used when Markowitz
1522 * product is to determine if off diagonals are placed symmetricly.
1523 * pOtherInRow (ElementPtr)
1524 * When there is only one other element in a row other than the diagonal,
1525 * pOtherInRow is used to point to it. Used when Markowitz product is
1526 * to determine if off diagonals are placed symmetricly.
1527 */
1528
1529 static ElementPtr
QuicklySearchDiagonal(Matrix,Step)1530 QuicklySearchDiagonal( Matrix, Step )
1531
1532 MatrixPtr Matrix;
1533 int Step;
1534 {
1535 register long MinMarkowitzProduct, *pMarkowitzProduct;
1536 register ElementPtr pDiag;
1537 int I;
1538 ElementPtr ChosenPivot, pOtherInRow, pOtherInCol;
1539 RealNumber Magnitude, LargestInCol, LargestOffDiagonal;
1540 RealNumber FindBiggestInColExclude();
1541
1542 /* Begin `QuicklySearchDiagonal'. */
1543 ChosenPivot = NULL;
1544 MinMarkowitzProduct = LARGEST_LONG_INTEGER;
1545 pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
1546 Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];
1547
1548 /* Assure that following while loop will always terminate. */
1549 Matrix->MarkowitzProd[Step-1] = -1;
1550
1551 /*
1552 * This is tricky. Am using a pointer in the inner while loop to
1553 * sequentially step through the MarkowitzProduct array. Search
1554 * terminates when the Markowitz product of zero placed at location
1555 * Step-1 is found. The row (and column) index on the diagonal is then
1556 * calculated by subtracting the pointer to the Markowitz product of
1557 * the first diagonal from the pointer to the Markowitz product of the
1558 * desired element. The outer for loop is infinite, broken by using
1559 * break.
1560 *
1561 * Search proceeds from the end (high row and column numbers) to the
1562 * beginning (low row and column numbers) so that rows and columns with
1563 * large Markowitz products will tend to be move to the bottom of the
1564 * matrix. However, choosing Diag[Step] is desirable because it would
1565 * require no row and column interchanges, so inspect it first by
1566 * putting its Markowitz product at the end of the MarkowitzProd
1567 * vector.
1568 */
1569
1570 for (;;) /* Endless for loop. */
1571 { while (*(--pMarkowitzProduct) >= MinMarkowitzProduct)
1572 { /* Just passing through. */
1573 }
1574
1575 I = pMarkowitzProduct - Matrix->MarkowitzProd;
1576
1577 /* Assure that I is valid; if I < Step, terminate search. */
1578 if (I < Step) break; /* Endless for loop */
1579 if (I > Matrix->Size) I = Step;
1580
1581 if ((pDiag = Matrix->Diag[I]) == NULL)
1582 continue; /* Endless for loop */
1583 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
1584 continue; /* Endless for loop */
1585
1586 if (*pMarkowitzProduct == 1)
1587 {
1588 /* Case where only one element exists in row and column other than diagonal. */
1589
1590 /* Find off-diagonal elements. */
1591 pOtherInRow = pDiag->NextInRow;
1592 pOtherInCol = pDiag->NextInCol;
1593 if (pOtherInRow == NULL AND pOtherInCol == NULL)
1594 { pOtherInRow = Matrix->FirstInRow[I];
1595 while(pOtherInRow != NULL)
1596 { if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I)
1597 break;
1598 pOtherInRow = pOtherInRow->NextInRow;
1599 }
1600 pOtherInCol = Matrix->FirstInCol[I];
1601 while(pOtherInCol != NULL)
1602 { if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I)
1603 break;
1604 pOtherInCol = pOtherInCol->NextInCol;
1605 }
1606 }
1607
1608 /* Accept diagonal as pivot if diagonal is larger than off-diagonals and the
1609 * off-diagonals are placed symmetricly. */
1610 if (pOtherInRow != NULL AND pOtherInCol != NULL)
1611 { if (pOtherInRow->Col == pOtherInCol->Row)
1612 { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
1613 ELEMENT_MAG(pOtherInCol));
1614 if (Magnitude >= LargestOffDiagonal)
1615 {
1616 /* Accept pivot, it is unlikely to contribute excess error. */
1617 return pDiag;
1618 }
1619 }
1620 }
1621 }
1622
1623 MinMarkowitzProduct = *pMarkowitzProduct;
1624 ChosenPivot = pDiag;
1625 } /* End of endless for loop. */
1626
1627 if (ChosenPivot != NULL)
1628 { LargestInCol = FindBiggestInColExclude( Matrix, ChosenPivot, Step );
1629 if( ELEMENT_MAG(ChosenPivot) <= Matrix->RelThreshold * LargestInCol )
1630 ChosenPivot = NULL;
1631 }
1632 return ChosenPivot;
1633 }
1634 #endif /* Not MODIFIED_MARKOWITZ */
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644 /*
1645 * SEARCH DIAGONAL FOR PIVOT WITH MODIFIED MARKOWITZ CRITERION
1646 *
1647 * Searches the diagonal looking for the best pivot. For a pivot to be
1648 * acceptable it must be larger than the pivot RelThreshold times the largest
1649 * element in its reduced column. Among the acceptable diagonals, the
1650 * one with the smallest MarkowitzProduct is sought. If a tie occurs
1651 * between elements of equal MarkowitzProduct, then the element with
1652 * the largest ratio between its magnitude and the largest element in its
1653 * column is used. The search will be terminated after a given number of
1654 * ties have occurred and the best (smallest ratio) of the tied element will
1655 * be used as the pivot. The number of ties that will trigger an early
1656 * termination is MinMarkowitzProduct * TIES_MULTIPLIER.
1657 *
1658 * >>> Returned:
1659 * A pointer to the diagonal element chosen to be pivot. If no diagonal is
1660 * acceptable, a NULL is returned.
1661 *
1662 * >>> Arguments:
1663 * Matrix <input> (MatrixPtr)
1664 * Pointer to matrix.
1665 * Step <input> (int)
1666 * Index of the diagonal currently being eliminated.
1667 *
1668 * >>> Local variables:
1669 * ChosenPivot (ElementPtr)
1670 * Pointer to the element that has been chosen to be the pivot.
1671 * Size (int)
1672 * Local version of size which is placed in a register to increase speed.
1673 * Magnitude (RealNumber)
1674 * Absolute value of diagonal element.
1675 * MinMarkowitzProduct (long)
1676 * Smallest Markowitz product found of those pivot candidates which are
1677 * acceptable.
1678 * NumberOfTies (int)
1679 * A count of the number of Markowitz ties that have occurred at current
1680 * MarkowitzProduct.
1681 * pDiag (ElementPtr)
1682 * Pointer to current diagonal element.
1683 * pMarkowitzProduct (long*)
1684 * Pointer that points into MarkowitzProduct array. It is used to quickly
1685 * access successive Markowitz products.
1686 * Ratio (RealNumber)
1687 * For the current pivot candidate, Ratio is the
1688 * Ratio of the largest element in its column to its magnitude.
1689 * RatioOfAccepted (RealNumber)
1690 * For the best pivot candidate found so far, RatioOfAccepted is the
1691 * Ratio of the largest element in its column to its magnitude.
1692 */
1693
1694 static ElementPtr
SearchDiagonal(Matrix,Step)1695 SearchDiagonal( Matrix, Step )
1696
1697 MatrixPtr Matrix;
1698 register int Step;
1699 {
1700 register int J;
1701 register long MinMarkowitzProduct, *pMarkowitzProduct;
1702 register int I;
1703 register ElementPtr pDiag;
1704 int NumberOfTies, Size = Matrix->Size;
1705 ElementPtr ChosenPivot;
1706 RealNumber Magnitude, Ratio, RatioOfAccepted, LargestInCol;
1707 RealNumber FindBiggestInColExclude();
1708
1709 /* Begin `SearchDiagonal'. */
1710 ChosenPivot = NULL;
1711 MinMarkowitzProduct = LARGEST_LONG_INTEGER;
1712 pMarkowitzProduct = &(Matrix->MarkowitzProd[Size+2]);
1713 Matrix->MarkowitzProd[Size+1] = Matrix->MarkowitzProd[Step];
1714
1715 /* Start search of diagonal. */
1716 for (J = Size+1; J > Step; J--)
1717 {
1718 if (*(--pMarkowitzProduct) > MinMarkowitzProduct)
1719 continue; /* for loop */
1720 if (J > Matrix->Size)
1721 I = Step;
1722 else
1723 I = J;
1724 if ((pDiag = Matrix->Diag[I]) == NULL)
1725 continue; /* for loop */
1726 if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
1727 continue; /* for loop */
1728
1729 /* Test to see if diagonal's magnitude is acceptable. */
1730 LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
1731 if (Magnitude <= Matrix->RelThreshold * LargestInCol)
1732 continue; /* for loop */
1733
1734 if (*pMarkowitzProduct < MinMarkowitzProduct)
1735 {
1736 /* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
1737 ChosenPivot = pDiag;
1738 MinMarkowitzProduct = *pMarkowitzProduct;
1739 RatioOfAccepted = LargestInCol / Magnitude;
1740 NumberOfTies = 0;
1741 }
1742 else
1743 {
1744 /* This case handles Markowitz ties. */
1745 NumberOfTies++;
1746 Ratio = LargestInCol / Magnitude;
1747 if (Ratio < RatioOfAccepted)
1748 { ChosenPivot = pDiag;
1749 RatioOfAccepted = Ratio;
1750 }
1751 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
1752 return ChosenPivot;
1753 }
1754 } /* End of for(Step) */
1755 return ChosenPivot;
1756 }
1757 #endif /* DIAGONAL_PIVOTING */
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768 /*
1769 * SEARCH ENTIRE MATRIX FOR BEST PIVOT
1770 *
1771 * Performs a search over the entire matrix looking for the acceptable
1772 * element with the lowest MarkowitzProduct. If there are several that
1773 * are tied for the smallest MarkowitzProduct, the tie is broken by using
1774 * the ratio of the magnitude of the element being considered to the largest
1775 * element in the same column. If no element is acceptable then the largest
1776 * element in the reduced submatrix is used as the pivot and the
1777 * matrix is declared to be spSMALL_PIVOT. If the largest element is
1778 * zero, the matrix is declared to be spSINGULAR.
1779 *
1780 * >>> Returned:
1781 * A pointer to the diagonal element chosen to be pivot. If no element is
1782 * found, then NULL is returned and the matrix is spSINGULAR.
1783 *
1784 * >>> Arguments:
1785 * Matrix <input> (MatrixPtr)
1786 * Pointer to matrix.
1787 * Step <input> (int)
1788 * Index of the diagonal currently being eliminated.
1789 *
1790 * >>> Local variables:
1791 * ChosenPivot (ElementPtr)
1792 * Pointer to the element that has been chosen to be the pivot.
1793 * LargestElementMag (RealNumber)
1794 * Magnitude of the largest element yet found in the reduced submatrix.
1795 * Size (int)
1796 * Local version of Size; placed in a register for speed.
1797 * Magnitude (RealNumber)
1798 * Absolute value of diagonal element.
1799 * MinMarkowitzProduct (long)
1800 * Smallest Markowitz product found of pivot candidates which are
1801 * acceptable.
1802 * NumberOfTies (int)
1803 * A count of the number of Markowitz ties that have occurred at current
1804 * MarkowitzProduct.
1805 * pElement (ElementPtr)
1806 * Pointer to current element.
1807 * pLargestElement (ElementPtr)
1808 * Pointer to the largest element yet found in the reduced submatrix.
1809 * Product (long)
1810 * Markowitz product for the current row and column.
1811 * Ratio (RealNumber)
1812 * For the current pivot candidate, Ratio is the
1813 * Ratio of the largest element in its column to its magnitude.
1814 * RatioOfAccepted (RealNumber)
1815 * For the best pivot candidate found so far, RatioOfAccepted is the
1816 * Ratio of the largest element in its column to its magnitude.
1817 *
1818 * >>> Possible errors:
1819 * spSINGULAR
1820 * spSMALL_PIVOT
1821 */
1822
1823 static ElementPtr
SearchEntireMatrix(Matrix,Step)1824 SearchEntireMatrix( Matrix, Step )
1825
1826 MatrixPtr Matrix;
1827 int Step;
1828 {
1829 register int I, Size = Matrix->Size;
1830 register ElementPtr pElement;
1831 int NumberOfTies;
1832 long Product, MinMarkowitzProduct;
1833 ElementPtr ChosenPivot, pLargestElement;
1834 RealNumber Magnitude, LargestElementMag, Ratio, RatioOfAccepted, LargestInCol;
1835 RealNumber FindLargestInCol();
1836
1837 /* Begin `SearchEntireMatrix'. */
1838 ChosenPivot = NULL;
1839 LargestElementMag = 0.0;
1840 MinMarkowitzProduct = LARGEST_LONG_INTEGER;
1841
1842 /* Start search of matrix on column by column basis. */
1843 for (I = Step; I <= Size; I++)
1844 { pElement = Matrix->FirstInCol[I];
1845
1846 while (pElement != NULL AND pElement->Row < Step)
1847 pElement = pElement->NextInCol;
1848
1849 if((LargestInCol = FindLargestInCol(pElement)) == 0.0)
1850 continue; /* for loop */
1851
1852 while (pElement != NULL)
1853 {
1854 /* Check to see if element is the largest encountered so far. If so, record
1855 its magnitude and address. */
1856 if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElementMag)
1857 { LargestElementMag = Magnitude;
1858 pLargestElement = pElement;
1859 }
1860 /* Calculate element's MarkowitzProduct. */
1861 Product = Matrix->MarkowitzRow[pElement->Row] *
1862 Matrix->MarkowitzCol[pElement->Col];
1863
1864 /* Test to see if element is acceptable as a pivot candidate. */
1865 if ((Product <= MinMarkowitzProduct) AND
1866 (Magnitude > Matrix->RelThreshold * LargestInCol) AND
1867 (Magnitude > Matrix->AbsThreshold))
1868 {
1869 /* Test to see if element has lowest MarkowitzProduct yet found, or whether it
1870 is tied with an element found earlier. */
1871 if (Product < MinMarkowitzProduct)
1872 {
1873 /* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
1874 ChosenPivot = pElement;
1875 MinMarkowitzProduct = Product;
1876 RatioOfAccepted = LargestInCol / Magnitude;
1877 NumberOfTies = 0;
1878 }
1879 else
1880 {
1881 /* This case handles Markowitz ties. */
1882 NumberOfTies++;
1883 Ratio = LargestInCol / Magnitude;
1884 if (Ratio < RatioOfAccepted)
1885 { ChosenPivot = pElement;
1886 RatioOfAccepted = Ratio;
1887 }
1888 if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
1889 return ChosenPivot;
1890 }
1891 }
1892 pElement = pElement->NextInCol;
1893 } /* End of while(pElement != NULL) */
1894 } /* End of for(Step) */
1895
1896 if (ChosenPivot != NULL) return ChosenPivot;
1897
1898 if (LargestElementMag == 0.0)
1899 { Matrix->Error = spSINGULAR;
1900 return NULL;
1901 }
1902
1903 Matrix->Error = spSMALL_PIVOT;
1904 return pLargestElement;
1905 }
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918 /*
1919 * DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
1920 *
1921 * This routine searches a column and returns the magnitude of the largest
1922 * element. This routine begins the search at the element pointed to by
1923 * pElement, the parameter.
1924 *
1925 * The search is conducted by starting at the element specified by a pointer,
1926 * which should be one below the diagonal, and moving down the column. On
1927 * the way down the column, the magnitudes of the elements are tested to see
1928 * if they are the largest yet found.
1929 *
1930 * >>> Returned:
1931 * The magnitude of the largest element in the column below and including
1932 * the one pointed to by the input parameter.
1933 *
1934 * >>> Arguments:
1935 * pElement <input> (ElementPtr)
1936 * The pointer to the first element to be tested. Also, used by the
1937 * routine to access all lower elements in the column.
1938 *
1939 * >>> Local variables:
1940 * Largest (RealNumber)
1941 * The magnitude of the largest element.
1942 * Magnitude (RealNumber)
1943 * The magnitude of the currently active element.
1944 */
1945
1946 static RealNumber
FindLargestInCol(pElement)1947 FindLargestInCol( pElement )
1948
1949 register ElementPtr pElement;
1950 {
1951 RealNumber Magnitude, Largest = 0.0;
1952
1953 /* Begin `FindLargestInCol'. */
1954 /* Search column for largest element beginning at Element. */
1955 while (pElement != NULL)
1956 { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
1957 Largest = Magnitude;
1958 pElement = pElement->NextInCol;
1959 }
1960
1961 return Largest;
1962 }
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973 /*
1974 * DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
1975 * EXCLUDING AN ELEMENT
1976 *
1977 * This routine searches a column and returns the magnitude of the largest
1978 * element. One given element is specifically excluded from the search.
1979 *
1980 * The search is conducted by starting at the first element in the column
1981 * and moving down the column until the active part of the matrix is entered,
1982 * i.e. the reduced submatrix. The rest of the column is then traversed
1983 * looking for the largest element.
1984 *
1985 * >>> Returned:
1986 * The magnitude of the largest element in the active portion of the column,
1987 * excluding the specified element, is returned.
1988 *
1989 * >>> Arguments:
1990 * Matrix <input> (MatrixPtr)
1991 * Pointer to the matrix.
1992 * pElement <input> (ElementPtr)
1993 * The pointer to the element that is to be excluded from search. Column
1994 * to be searched is one that contains this element. Also used to
1995 * access the elements in the column.
1996 * Step <input> (int)
1997 * Index of the diagonal currently being eliminated. Indicates where
1998 * the active part of the matrix begins.
1999 *
2000 * >>> Local variables:
2001 * Col (int)
2002 * The number of the column to be searched. Also the column number of
2003 * the element to be avoided in the search.
2004 * Largest (RealNumber)
2005 * The magnitude of the largest element.
2006 * Magnitude (RealNumber)
2007 * The magnitude of the currently active element.
2008 * Row (int)
2009 * The row number of element to be excluded from the search.
2010 */
2011
2012 static RealNumber
FindBiggestInColExclude(Matrix,pElement,Step)2013 FindBiggestInColExclude( Matrix, pElement, Step )
2014
2015 MatrixPtr Matrix;
2016 register ElementPtr pElement;
2017 register int Step;
2018 {
2019 register int Row;
2020 int Col;
2021 RealNumber Largest, Magnitude;
2022
2023 /* Begin `FindBiggestInColExclude'. */
2024 Row = pElement->Row;
2025 Col = pElement->Col;
2026 pElement = Matrix->FirstInCol[Col];
2027
2028 /* Travel down column until reduced submatrix is entered. */
2029 while ((pElement != NULL) AND (pElement->Row < Step))
2030 pElement = pElement->NextInCol;
2031
2032 /* Initialize the variable Largest. */
2033 if (pElement->Row != Row)
2034 Largest = ELEMENT_MAG(pElement);
2035 else
2036 Largest = 0.0;
2037
2038 /* Search rest of column for largest element, avoiding excluded element. */
2039 while ((pElement = pElement->NextInCol) != NULL)
2040 { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
2041 { if (pElement->Row != Row)
2042 Largest = Magnitude;
2043 }
2044 }
2045
2046 return Largest;
2047 }
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058 /*
2059 * EXCHANGE ROWS AND COLUMNS
2060 *
2061 * Exchanges two rows and two columns so that the selected pivot is moved to
2062 * the upper left corner of the remaining submatrix.
2063 *
2064 * >>> Arguments:
2065 * Matrix <input> (MatrixPtr)
2066 * Pointer to the matrix.
2067 * pPivot <input> (ElementPtr)
2068 * Pointer to the current pivot.
2069 * Step <input> (int)
2070 * Index of the diagonal currently being eliminated.
2071 *
2072 * >>> Local variables:
2073 * Col (int)
2074 * Column where the pivot was found.
2075 * Row (int)
2076 * Row where the pivot was found.
2077 * OldMarkowitzProd_Col (long)
2078 * Markowitz product associated with the diagonal element in the row
2079 * the pivot was found in.
2080 * OldMarkowitzProd_Row (long)
2081 * Markowitz product associated with the diagonal element in the column
2082 * the pivot was found in.
2083 * OldMarkowitzProd_Step (long)
2084 * Markowitz product associated with the diagonal element that is being
2085 * moved so that the pivot can be placed in the upper left-hand corner
2086 * of the reduced submatrix.
2087 */
2088
2089 static
ExchangeRowsAndCols(Matrix,pPivot,Step)2090 ExchangeRowsAndCols( Matrix, pPivot, Step )
2091
2092 MatrixPtr Matrix;
2093 ElementPtr pPivot;
2094 register int Step;
2095 {
2096 register int Row, Col;
2097 long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col;
2098 ElementPtr spcFindElementInCol();
2099
2100 /* Begin `ExchangeRowsAndCols'. */
2101 Row = pPivot->Row;
2102 Col = pPivot->Col;
2103 Matrix->PivotsOriginalRow = Row;
2104 Matrix->PivotsOriginalCol = Col;
2105
2106 if ((Row == Step) AND (Col == Step)) return 0;
2107
2108 /* Exchange rows and columns. */
2109 if (Row == Col)
2110 { spcRowExchange( Matrix, Step, Row );
2111 spcColExchange( Matrix, Step, Col );
2112 SWAP( long, Matrix->MarkowitzProd[Step], Matrix->MarkowitzProd[Row] );
2113 SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Step] );
2114 }
2115 else
2116 {
2117
2118 /* Initialize variables that hold old Markowitz products. */
2119 OldMarkowitzProd_Step = Matrix->MarkowitzProd[Step];
2120 OldMarkowitzProd_Row = Matrix->MarkowitzProd[Row];
2121 OldMarkowitzProd_Col = Matrix->MarkowitzProd[Col];
2122
2123 /* Exchange rows. */
2124 if (Row != Step)
2125 { spcRowExchange( Matrix, Step, Row );
2126 Matrix->NumberOfInterchangesIsOdd =
2127 NOT Matrix->NumberOfInterchangesIsOdd;
2128 Matrix->MarkowitzProd[Row] = Matrix->MarkowitzRow[Row] *
2129 Matrix->MarkowitzCol[Row];
2130
2131 /* Update singleton count. */
2132 if ((Matrix->MarkowitzProd[Row]==0) != (OldMarkowitzProd_Row==0))
2133 { if (OldMarkowitzProd_Row == 0)
2134 Matrix->Singletons--;
2135 else
2136 Matrix->Singletons++;
2137 }
2138 }
2139
2140 /* Exchange columns. */
2141 if (Col != Step)
2142 { spcColExchange( Matrix, Step, Col );
2143 Matrix->NumberOfInterchangesIsOdd =
2144 NOT Matrix->NumberOfInterchangesIsOdd;
2145 Matrix->MarkowitzProd[Col] = Matrix->MarkowitzCol[Col] *
2146 Matrix->MarkowitzRow[Col];
2147
2148 /* Update singleton count. */
2149 if ((Matrix->MarkowitzProd[Col]==0) != (OldMarkowitzProd_Col==0))
2150 { if (OldMarkowitzProd_Col == 0)
2151 Matrix->Singletons--;
2152 else
2153 Matrix->Singletons++;
2154 }
2155
2156 Matrix->Diag[Col] = spcFindElementInCol( Matrix,
2157 Matrix->FirstInCol+Col,
2158 Col, Col, NO );
2159 }
2160 if (Row != Step)
2161 { Matrix->Diag[Row] = spcFindElementInCol( Matrix,
2162 Matrix->FirstInCol+Row,
2163 Row, Row, NO );
2164 }
2165 Matrix->Diag[Step] = spcFindElementInCol( Matrix,
2166 Matrix->FirstInCol+Step,
2167 Step, Step, NO );
2168
2169 /* Update singleton count. */
2170 Matrix->MarkowitzProd[Step] = Matrix->MarkowitzCol[Step] *
2171 Matrix->MarkowitzRow[Step];
2172 if ((Matrix->MarkowitzProd[Step]==0) != (OldMarkowitzProd_Step==0))
2173 { if (OldMarkowitzProd_Step == 0)
2174 Matrix->Singletons--;
2175 else
2176 Matrix->Singletons++;
2177 }
2178 }
2179 return 0;
2180 }
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190 /*
2191 * EXCHANGE ROWS
2192 *
2193 * Performs all required operations to exchange two rows. Those operations
2194 * include: swap FirstInRow pointers, fixing up the NextInCol pointers,
2195 * swapping row indexes in MatrixElements, and swapping Markowitz row
2196 * counts.
2197 *
2198 * >>> Arguments:
2199 * Matrix <input> (MatrixPtr)
2200 * Pointer to the matrix.
2201 * Row1 <input> (int)
2202 * Row index of one of the rows, becomes the smallest index.
2203 * Row2 <input> (int)
2204 * Row index of the other row, becomes the largest index.
2205 *
2206 * Local variables:
2207 * Column (int)
2208 * Column in which row elements are currently being exchanged.
2209 * Row1Ptr (ElementPtr)
2210 * Pointer to an element in Row1.
2211 * Row2Ptr (ElementPtr)
2212 * Pointer to an element in Row2.
2213 * Element1 (ElementPtr)
2214 * Pointer to the element in Row1 to be exchanged.
2215 * Element2 (ElementPtr)
2216 * Pointer to the element in Row2 to be exchanged.
2217 */
2218
spcRowExchange(Matrix,Row1,Row2)2219 spcRowExchange( Matrix, Row1, Row2 )
2220
2221 MatrixPtr Matrix;
2222 int Row1, Row2;
2223 {
2224 register ElementPtr Row1Ptr, Row2Ptr;
2225 int Column;
2226 ElementPtr Element1, Element2;
2227
2228 /* Begin `spcRowExchange'. */
2229 if (Row1 > Row2) SWAP(int, Row1, Row2);
2230
2231 Row1Ptr = Matrix->FirstInRow[Row1];
2232 Row2Ptr = Matrix->FirstInRow[Row2];
2233 while (Row1Ptr != NULL OR Row2Ptr != NULL)
2234 {
2235 /* Exchange elements in rows while traveling from left to right. */
2236 if (Row1Ptr == NULL)
2237 { Column = Row2Ptr->Col;
2238 Element1 = NULL;
2239 Element2 = Row2Ptr;
2240 Row2Ptr = Row2Ptr->NextInRow;
2241 }
2242 else if (Row2Ptr == NULL)
2243 { Column = Row1Ptr->Col;
2244 Element1 = Row1Ptr;
2245 Element2 = NULL;
2246 Row1Ptr = Row1Ptr->NextInRow;
2247 }
2248 else if (Row1Ptr->Col < Row2Ptr->Col)
2249 { Column = Row1Ptr->Col;
2250 Element1 = Row1Ptr;
2251 Element2 = NULL;
2252 Row1Ptr = Row1Ptr->NextInRow;
2253 }
2254 else if (Row1Ptr->Col > Row2Ptr->Col)
2255 { Column = Row2Ptr->Col;
2256 Element1 = NULL;
2257 Element2 = Row2Ptr;
2258 Row2Ptr = Row2Ptr->NextInRow;
2259 }
2260 else /* Row1Ptr->Col == Row2Ptr->Col */
2261 { Column = Row1Ptr->Col;
2262 Element1 = Row1Ptr;
2263 Element2 = Row2Ptr;
2264 Row1Ptr = Row1Ptr->NextInRow;
2265 Row2Ptr = Row2Ptr->NextInRow;
2266 }
2267
2268 ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column);
2269 } /* end of while(Row1Ptr != NULL OR Row2Ptr != NULL) */
2270
2271 if (Matrix->InternalVectorsAllocated)
2272 SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]);
2273 SWAP( ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]);
2274 SWAP( int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]);
2275 #if TRANSLATE
2276 Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row1] ] = Row1;
2277 Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row2] ] = Row2;
2278 #endif
2279
2280 return 0;
2281 }
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291 /*
2292 * EXCHANGE COLUMNS
2293 *
2294 * Performs all required operations to exchange two columns. Those operations
2295 * include: swap FirstInCol pointers, fixing up the NextInRow pointers,
2296 * swapping column indexes in MatrixElements, and swapping Markowitz
2297 * column counts.
2298 *
2299 * >>> Arguments:
2300 * Matrix <input> (MatrixPtr)
2301 * Pointer to the matrix.
2302 * Col1 <input> (int)
2303 * Column index of one of the columns, becomes the smallest index.
2304 * Col2 <input> (int)
2305 * Column index of the other column, becomes the largest index
2306 *
2307 * Local variables:
2308 * Row (int)
2309 * Row in which column elements are currently being exchanged.
2310 * Col1Ptr (ElementPtr)
2311 * Pointer to an element in Col1.
2312 * Col2Ptr (ElementPtr)
2313 * Pointer to an element in Col2.
2314 * Element1 (ElementPtr)
2315 * Pointer to the element in Col1 to be exchanged.
2316 * Element2 (ElementPtr)
2317 * Pointer to the element in Col2 to be exchanged.
2318 */
2319
spcColExchange(Matrix,Col1,Col2)2320 spcColExchange( Matrix, Col1, Col2 )
2321
2322 MatrixPtr Matrix;
2323 int Col1, Col2;
2324 {
2325 register ElementPtr Col1Ptr, Col2Ptr;
2326 int Row;
2327 ElementPtr Element1, Element2;
2328
2329 /* Begin `spcColExchange'. */
2330 if (Col1 > Col2) SWAP(int, Col1, Col2);
2331
2332 Col1Ptr = Matrix->FirstInCol[Col1];
2333 Col2Ptr = Matrix->FirstInCol[Col2];
2334 while (Col1Ptr != NULL OR Col2Ptr != NULL)
2335 {
2336 /* Exchange elements in rows while traveling from top to bottom. */
2337 if (Col1Ptr == NULL)
2338 { Row = Col2Ptr->Row;
2339 Element1 = NULL;
2340 Element2 = Col2Ptr;
2341 Col2Ptr = Col2Ptr->NextInCol;
2342 }
2343 else if (Col2Ptr == NULL)
2344 { Row = Col1Ptr->Row;
2345 Element1 = Col1Ptr;
2346 Element2 = NULL;
2347 Col1Ptr = Col1Ptr->NextInCol;
2348 }
2349 else if (Col1Ptr->Row < Col2Ptr->Row)
2350 { Row = Col1Ptr->Row;
2351 Element1 = Col1Ptr;
2352 Element2 = NULL;
2353 Col1Ptr = Col1Ptr->NextInCol;
2354 }
2355 else if (Col1Ptr->Row > Col2Ptr->Row)
2356 { Row = Col2Ptr->Row;
2357 Element1 = NULL;
2358 Element2 = Col2Ptr;
2359 Col2Ptr = Col2Ptr->NextInCol;
2360 }
2361 else /* Col1Ptr->Row == Col2Ptr->Row */
2362 { Row = Col1Ptr->Row;
2363 Element1 = Col1Ptr;
2364 Element2 = Col2Ptr;
2365 Col1Ptr = Col1Ptr->NextInCol;
2366 Col2Ptr = Col2Ptr->NextInCol;
2367 }
2368
2369 ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row);
2370 } /* end of while(Col1Ptr != NULL OR Col2Ptr != NULL) */
2371
2372 if (Matrix->InternalVectorsAllocated)
2373 SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]);
2374 SWAP( ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
2375 SWAP( int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
2376 #if TRANSLATE
2377 Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col1] ] = Col1;
2378 Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col2] ] = Col2;
2379 #endif
2380
2381 return 0;
2382 }
2383
2384
2385
2386
2387
2388
2389
2390 /*
2391 * EXCHANGE TWO ELEMENTS IN A COLUMN
2392 *
2393 * Performs all required operations to exchange two elements in a column.
2394 * Those operations are: restring NextInCol pointers and swapping row indexes
2395 * in the MatrixElements.
2396 *
2397 * >>> Arguments:
2398 * Matrix <input> (MatrixPtr)
2399 * Pointer to the matrix.
2400 * Row1 <input> (int)
2401 * Row of top element to be exchanged.
2402 * Element1 <input> (ElementPtr)
2403 * Pointer to top element to be exchanged.
2404 * Row2 <input> (int)
2405 * Row of bottom element to be exchanged.
2406 * Element2 <input> (ElementPtr)
2407 * Pointer to bottom element to be exchanged.
2408 * Column <input> (int)
2409 * Column that exchange is to take place in.
2410 *
2411 * >>> Local variables:
2412 * ElementAboveRow1 (ElementPtr *)
2413 * Location of pointer which points to the element above Element1. This
2414 * pointer is modified so that it points to correct element on exit.
2415 * ElementAboveRow2 (ElementPtr *)
2416 * Location of pointer which points to the element above Element2. This
2417 * pointer is modified so that it points to correct element on exit.
2418 * ElementBelowRow1 (ElementPtr)
2419 * Pointer to element below Element1.
2420 * ElementBelowRow2 (ElementPtr)
2421 * Pointer to element below Element2.
2422 * pElement (ElementPtr)
2423 * Pointer used to traverse the column.
2424 */
2425
2426 static
ExchangeColElements(Matrix,Row1,Element1,Row2,Element2,Column)2427 ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column )
2428
2429 MatrixPtr Matrix;
2430 register ElementPtr Element1, Element2;
2431 int Row1, Row2, Column;
2432 {
2433 ElementPtr *ElementAboveRow1, *ElementAboveRow2;
2434 ElementPtr ElementBelowRow1, ElementBelowRow2;
2435 register ElementPtr pElement;
2436
2437 /* Begin `ExchangeColElements'. */
2438 /* Search to find the ElementAboveRow1. */
2439 ElementAboveRow1 = &(Matrix->FirstInCol[Column]);
2440 pElement = *ElementAboveRow1;
2441 while (pElement->Row < Row1)
2442 { ElementAboveRow1 = &(pElement->NextInCol);
2443 pElement = *ElementAboveRow1;
2444 }
2445 if (Element1 != NULL)
2446 { ElementBelowRow1 = Element1->NextInCol;
2447 if (Element2 == NULL)
2448 {
2449 /* Element2 does not exist, move Element1 down to Row2. */
2450 if ( ElementBelowRow1 != NULL AND ElementBelowRow1->Row < Row2 )
2451 {
2452 /* Element1 must be removed from linked list and moved. */
2453 *ElementAboveRow1 = ElementBelowRow1;
2454
2455 /* Search column for Row2. */
2456 pElement = ElementBelowRow1;
2457 do
2458 { ElementAboveRow2 = &(pElement->NextInCol);
2459 pElement = *ElementAboveRow2;
2460 } while (pElement != NULL AND pElement->Row < Row2);
2461
2462 /* Place Element1 in Row2. */
2463 *ElementAboveRow2 = Element1;
2464 Element1->NextInCol = pElement;
2465 *ElementAboveRow1 =ElementBelowRow1;
2466 }
2467 Element1->Row = Row2;
2468 }
2469 else
2470 {
2471 /* Element2 does exist, and the two elements must be exchanged. */
2472 if ( ElementBelowRow1->Row == Row2)
2473 {
2474 /* Element2 is just below Element1, exchange them. */
2475 Element1->NextInCol = Element2->NextInCol;
2476 Element2->NextInCol = Element1;
2477 *ElementAboveRow1 = Element2;
2478 }
2479 else
2480 {
2481 /* Element2 is not just below Element1 and must be searched for. */
2482 pElement = ElementBelowRow1;
2483 do
2484 { ElementAboveRow2 = &(pElement->NextInCol);
2485 pElement = *ElementAboveRow2;
2486 } while (pElement->Row < Row2);
2487
2488 ElementBelowRow2 = Element2->NextInCol;
2489
2490 /* Switch Element1 and Element2. */
2491 *ElementAboveRow1 = Element2;
2492 Element2->NextInCol = ElementBelowRow1;
2493 *ElementAboveRow2 = Element1;
2494 Element1->NextInCol = ElementBelowRow2;
2495 }
2496 Element1->Row = Row2;
2497 Element2->Row = Row1;
2498 }
2499 }
2500 else
2501 {
2502 /* Element1 does not exist. */
2503 ElementBelowRow1 = pElement;
2504
2505 /* Find Element2. */
2506 if (ElementBelowRow1->Row != Row2)
2507 { do
2508 { ElementAboveRow2 = &(pElement->NextInCol);
2509 pElement = *ElementAboveRow2;
2510 } while (pElement->Row < Row2);
2511
2512 ElementBelowRow2 = Element2->NextInCol;
2513
2514 /* Move Element2 to Row1. */
2515 *ElementAboveRow2 = Element2->NextInCol;
2516 *ElementAboveRow1 = Element2;
2517 Element2->NextInCol = ElementBelowRow1;
2518 }
2519 Element2->Row = Row1;
2520 }
2521 return 0;
2522 }
2523
2524
2525
2526
2527
2528
2529
2530 /*
2531 * EXCHANGE TWO ELEMENTS IN A ROW
2532 *
2533 * Performs all required operations to exchange two elements in a row.
2534 * Those operations are: restring NextInRow pointers and swapping column
2535 * indexes in the MatrixElements.
2536 *
2537 * >>> Arguments:
2538 * Matrix <input> (MatrixPtr)
2539 * Pointer to the matrix.
2540 * Col1 <input> (int)
2541 * Col of left-most element to be exchanged.
2542 * Element1 <input> (ElementPtr)
2543 * Pointer to left-most element to be exchanged.
2544 * Col2 <input> (int)
2545 * Col of right-most element to be exchanged.
2546 * Element2 <input> (ElementPtr)
2547 * Pointer to right-most element to be exchanged.
2548 * Row <input> (int)
2549 * Row that exchange is to take place in.
2550 *
2551 * >>> Local variables:
2552 * ElementLeftOfCol1 (ElementPtr *)
2553 * Location of pointer which points to the element to the left of
2554 * Element1. This pointer is modified so that it points to correct
2555 * element on exit.
2556 * ElementLeftOfCol2 (ElementPtr *)
2557 * Location of pointer which points to the element to the left of
2558 * Element2. This pointer is modified so that it points to correct
2559 * element on exit.
2560 * ElementRightOfCol1 (ElementPtr)
2561 * Pointer to element right of Element1.
2562 * ElementRightOfCol2 (ElementPtr)
2563 * Pointer to element right of Element2.
2564 * pElement (ElementPtr)
2565 * Pointer used to traverse the row.
2566 */
2567
2568 static
ExchangeRowElements(Matrix,Col1,Element1,Col2,Element2,Row)2569 ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row )
2570
2571 MatrixPtr Matrix;
2572 int Col1, Col2, Row;
2573 register ElementPtr Element1, Element2;
2574 {
2575 ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2;
2576 ElementPtr ElementRightOfCol1, ElementRightOfCol2;
2577 register ElementPtr pElement;
2578
2579 /* Begin `ExchangeRowElements'. */
2580 /* Search to find the ElementLeftOfCol1. */
2581 ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]);
2582 pElement = *ElementLeftOfCol1;
2583 while (pElement->Col < Col1)
2584 { ElementLeftOfCol1 = &(pElement->NextInRow);
2585 pElement = *ElementLeftOfCol1;
2586 }
2587 if (Element1 != NULL)
2588 { ElementRightOfCol1 = Element1->NextInRow;
2589 if (Element2 == NULL)
2590 {
2591 /* Element2 does not exist, move Element1 to right to Col2. */
2592 if ( ElementRightOfCol1 != NULL AND ElementRightOfCol1->Col < Col2 )
2593 {
2594 /* Element1 must be removed from linked list and moved. */
2595 *ElementLeftOfCol1 = ElementRightOfCol1;
2596
2597 /* Search Row for Col2. */
2598 pElement = ElementRightOfCol1;
2599 do
2600 { ElementLeftOfCol2 = &(pElement->NextInRow);
2601 pElement = *ElementLeftOfCol2;
2602 } while (pElement != NULL AND pElement->Col < Col2);
2603
2604 /* Place Element1 in Col2. */
2605 *ElementLeftOfCol2 = Element1;
2606 Element1->NextInRow = pElement;
2607 *ElementLeftOfCol1 =ElementRightOfCol1;
2608 }
2609 Element1->Col = Col2;
2610 }
2611 else
2612 {
2613 /* Element2 does exist, and the two elements must be exchanged. */
2614 if ( ElementRightOfCol1->Col == Col2)
2615 {
2616 /* Element2 is just right of Element1, exchange them. */
2617 Element1->NextInRow = Element2->NextInRow;
2618 Element2->NextInRow = Element1;
2619 *ElementLeftOfCol1 = Element2;
2620 }
2621 else
2622 {
2623 /* Element2 is not just right of Element1 and must be searched for. */
2624 pElement = ElementRightOfCol1;
2625 do
2626 { ElementLeftOfCol2 = &(pElement->NextInRow);
2627 pElement = *ElementLeftOfCol2;
2628 } while (pElement->Col < Col2);
2629
2630 ElementRightOfCol2 = Element2->NextInRow;
2631
2632 /* Switch Element1 and Element2. */
2633 *ElementLeftOfCol1 = Element2;
2634 Element2->NextInRow = ElementRightOfCol1;
2635 *ElementLeftOfCol2 = Element1;
2636 Element1->NextInRow = ElementRightOfCol2;
2637 }
2638 Element1->Col = Col2;
2639 Element2->Col = Col1;
2640 }
2641 }
2642 else
2643 {
2644 /* Element1 does not exist. */
2645 ElementRightOfCol1 = pElement;
2646
2647 /* Find Element2. */
2648 if (ElementRightOfCol1->Col != Col2)
2649 { do
2650 { ElementLeftOfCol2 = &(pElement->NextInRow);
2651 pElement = *ElementLeftOfCol2;
2652 } while (pElement->Col < Col2);
2653
2654 ElementRightOfCol2 = Element2->NextInRow;
2655
2656 /* Move Element2 to Col1. */
2657 *ElementLeftOfCol2 = Element2->NextInRow;
2658 *ElementLeftOfCol1 = Element2;
2659 Element2->NextInRow = ElementRightOfCol1;
2660 }
2661 Element2->Col = Col1;
2662 }
2663 return 0;
2664 }
2665
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676 /*
2677 * PERFORM ROW AND COLUMN ELIMINATION ON REAL MATRIX
2678 *
2679 * Eliminates a single row and column of the matrix and leaves single row of
2680 * the upper triangular matrix and a single column of the lower triangular
2681 * matrix in its wake. Uses Gauss's method.
2682 *
2683 * >>> Argument:
2684 * Matrix <input> (MatrixPtr)
2685 * Pointer to the matrix.
2686 * pPivot <input> (ElementPtr)
2687 * Pointer to the current pivot.
2688 *
2689 * >>> Local variables:
2690 * pLower (ElementPtr)
2691 * Points to matrix element in lower triangular column.
2692 * pSub (ElementPtr)
2693 * Points to elements in the reduced submatrix.
2694 * Row (int)
2695 * Row index.
2696 * pUpper (ElementPtr)
2697 * Points to matrix element in upper triangular row.
2698 *
2699 * >>> Possible errors:
2700 * spNO_MEMORY
2701 */
2702
2703 static
RealRowColElimination(Matrix,pPivot)2704 RealRowColElimination( Matrix, pPivot )
2705
2706 MatrixPtr Matrix;
2707 register ElementPtr pPivot;
2708 {
2709 #if REAL
2710 register ElementPtr pSub;
2711 register int Row;
2712 register ElementPtr pLower, pUpper;
2713 extern ElementPtr CreateFillin();
2714
2715 /* Begin `RealRowColElimination'. */
2716
2717 /* Test for zero pivot. */
2718 if (ABS(pPivot->Real) == 0.0)
2719 { (void)MatrixIsSingular( Matrix, pPivot->Row );
2720 return 0;
2721 }
2722 pPivot->Real = 1.0 / pPivot->Real;
2723
2724 pUpper = pPivot->NextInRow;
2725 while (pUpper != NULL)
2726 {
2727 /* Calculate upper triangular element. */
2728 pUpper->Real *= pPivot->Real;
2729
2730 pSub = pUpper->NextInCol;
2731 pLower = pPivot->NextInCol;
2732 while (pLower != NULL)
2733 { Row = pLower->Row;
2734
2735 /* Find element in row that lines up with current lower triangular element. */
2736 while (pSub != NULL AND pSub->Row < Row)
2737 pSub = pSub->NextInCol;
2738
2739 /* Test to see if desired element was not found, if not, create fill-in. */
2740 if (pSub == NULL OR pSub->Row > Row)
2741 { pSub = CreateFillin( Matrix, Row, pUpper->Col );
2742 if (pSub == NULL)
2743 { Matrix->Error = spNO_MEMORY;
2744 return 0;
2745 }
2746 }
2747 pSub->Real -= pUpper->Real * pLower->Real;
2748 pSub = pSub->NextInCol;
2749 pLower = pLower->NextInCol;
2750 }
2751 pUpper = pUpper->NextInRow;
2752 }
2753 return 0;
2754 #endif /* REAL */
2755 }
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765 /*
2766 * PERFORM ROW AND COLUMN ELIMINATION ON COMPLEX MATRIX
2767 *
2768 * Eliminates a single row and column of the matrix and leaves single row of
2769 * the upper triangular matrix and a single column of the lower triangular
2770 * matrix in its wake. Uses Gauss's method.
2771 *
2772 * >>> Argument:
2773 * Matrix <input> (MatrixPtr)
2774 * Pointer to the matrix.
2775 * pPivot <input> (ElementPtr)
2776 * Pointer to the current pivot.
2777 *
2778 * >>> Local variables:
2779 * pLower (ElementPtr)
2780 * Points to matrix element in lower triangular column.
2781 * pSub (ElementPtr)
2782 * Points to elements in the reduced submatrix.
2783 * Row (int)
2784 * Row index.
2785 * pUpper (ElementPtr)
2786 * Points to matrix element in upper triangular row.
2787 *
2788 * Possible errors:
2789 * spNO_MEMORY
2790 */
2791
2792 static
ComplexRowColElimination(Matrix,pPivot)2793 ComplexRowColElimination( Matrix, pPivot )
2794
2795 MatrixPtr Matrix;
2796 register ElementPtr pPivot;
2797 {
2798 #if spCOMPLEX
2799 register ElementPtr pSub;
2800 register int Row;
2801 register ElementPtr pLower, pUpper;
2802 ElementPtr CreateFillin();
2803
2804 /* Begin `ComplexRowColElimination'. */
2805
2806 /* Test for zero pivot. */
2807 if (ELEMENT_MAG(pPivot) == 0.0)
2808 { (void)MatrixIsSingular( Matrix, pPivot->Row );
2809 return 0;
2810 }
2811 CMPLX_RECIPROCAL(*pPivot, *pPivot);
2812
2813 pUpper = pPivot->NextInRow;
2814 while (pUpper != NULL)
2815 {
2816 /* Calculate upper triangular element. */
2817 /* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */
2818 CMPLX_MULT_ASSIGN(*pUpper, *pPivot);
2819
2820 pSub = pUpper->NextInCol;
2821 pLower = pPivot->NextInCol;
2822 while (pLower != NULL)
2823 { Row = pLower->Row;
2824
2825 /* Find element in row that lines up with current lower triangular element. */
2826 while (pSub != NULL AND pSub->Row < Row)
2827 pSub = pSub->NextInCol;
2828
2829 /* Test to see if desired element was not found, if not, create fill-in. */
2830 if (pSub == NULL OR pSub->Row > Row)
2831 { pSub = CreateFillin( Matrix, Row, pUpper->Col );
2832 if (pSub == NULL)
2833 { Matrix->Error = spNO_MEMORY;
2834 return 0;
2835 }
2836 }
2837
2838 /* Cmplx expr: pElement -= *pUpper * pLower. */
2839 CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower);
2840 pSub = pSub->NextInCol;
2841 pLower = pLower->NextInCol;
2842 }
2843 pUpper = pUpper->NextInRow;
2844 }
2845 return 0;
2846 #endif /* spCOMPLEX */
2847 }
2848
2849
2850
2851
2852
2853 /*
2854 * UPDATE MARKOWITZ NUMBERS
2855 *
2856 * Updates the Markowitz numbers after a row and column have been eliminated.
2857 * Also updates singleton count.
2858 *
2859 * >>> Argument:
2860 * Matrix <input> (MatrixPtr)
2861 * Pointer to the matrix.
2862 * pPivot <input> (ElementPtr)
2863 * Pointer to the current pivot.
2864 *
2865 * >>> Local variables:
2866 * Row (int)
2867 * Row index.
2868 * Col (int)
2869 * Column index.
2870 * ColPtr (ElementPtr)
2871 * Points to matrix element in upper triangular column.
2872 * RowPtr (ElementPtr)
2873 * Points to matrix element in lower triangular row.
2874 */
2875
2876 static
UpdateMarkowitzNumbers(Matrix,pPivot)2877 UpdateMarkowitzNumbers( Matrix, pPivot )
2878
2879 MatrixPtr Matrix;
2880 ElementPtr pPivot;
2881 {
2882 register int Row, Col;
2883 register ElementPtr ColPtr, RowPtr;
2884 register int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
2885 double Product;
2886
2887 /* Begin `UpdateMarkowitzNumbers'. */
2888
2889 /* Update Markowitz numbers. */
2890 for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol)
2891 { Row = ColPtr->Row;
2892 --MarkoRow[Row];
2893
2894 /* Form Markowitz product while being cautious of overflows. */
2895 if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER AND MarkoCol[Row] != 0) OR
2896 (MarkoCol[Row] > LARGEST_SHORT_INTEGER AND MarkoRow[Row] != 0))
2897 { Product = MarkoCol[Row] * MarkoRow[Row];
2898 if (Product >= LARGEST_LONG_INTEGER)
2899 Matrix->MarkowitzProd[Row] = LARGEST_LONG_INTEGER;
2900 else
2901 Matrix->MarkowitzProd[Row] = Product;
2902 }
2903 else Matrix->MarkowitzProd[Row] = MarkoRow[Row] * MarkoCol[Row];
2904 if (MarkoRow[Row] == 0)
2905 Matrix->Singletons++;
2906 }
2907
2908 for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow)
2909 { Col = RowPtr->Col;
2910 --MarkoCol[Col];
2911
2912 /* Form Markowitz product while being cautious of overflows. */
2913 if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER AND MarkoCol[Col] != 0) OR
2914 (MarkoCol[Col] > LARGEST_SHORT_INTEGER AND MarkoRow[Col] != 0))
2915 { Product = MarkoCol[Col] * MarkoRow[Col];
2916 if (Product >= LARGEST_LONG_INTEGER)
2917 Matrix->MarkowitzProd[Col] = LARGEST_LONG_INTEGER;
2918 else
2919 Matrix->MarkowitzProd[Col] = Product;
2920 }
2921 else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col];
2922 if ((MarkoCol[Col] == 0) AND (MarkoRow[Col] != 0))
2923 Matrix->Singletons++;
2924 }
2925 return 0;
2926 }
2927
2928
2929
2930
2931
2932
2933
2934
2935 /*
2936 * CREATE FILL-IN
2937 *
2938 * This routine is used to create fill-ins and splice them into the
2939 * matrix.
2940 *
2941 * >>> Returns:
2942 * Pointer to fill-in.
2943 *
2944 * >>> Arguments:
2945 * Matrix <input> (MatrixPtr)
2946 * Pointer to the matrix.
2947 * Col <input> (int)
2948 * Column index for element.
2949 * Row <input> (int)
2950 * Row index for element.
2951 *
2952 * >>> Local variables:
2953 * pElement (ElementPtr)
2954 * Pointer to an element in the matrix.
2955 * ppElementAbove (ElementPtr *)
2956 * This contains the address of the pointer to the element just above the
2957 * one being created. It is used to speed the search and it is updated
2958 * with address of the created element.
2959 *
2960 * >>> Possible errors:
2961 * spNO_MEMORY
2962 */
2963
2964 static ElementPtr
CreateFillin(Matrix,Row,Col)2965 CreateFillin( Matrix, Row, Col )
2966
2967 MatrixPtr Matrix;
2968 register int Row;
2969 int Col;
2970 {
2971 register ElementPtr pElement, *ppElementAbove;
2972 ElementPtr spcCreateElement();
2973
2974 /* Begin `CreateFillin'. */
2975
2976 /* Find Element above fill-in. */
2977 ppElementAbove = &Matrix->FirstInCol[Col];
2978 pElement = *ppElementAbove;
2979 while (pElement != NULL)
2980 { if (pElement->Row < Row)
2981 { ppElementAbove = &pElement->NextInCol;
2982 pElement = *ppElementAbove;
2983 }
2984 else break; /* while loop */
2985 }
2986
2987 /* End of search, create the element. */
2988 pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, YES );
2989
2990 /* Update Markowitz counts and products. */
2991 Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] *
2992 Matrix->MarkowitzCol[Row];
2993 if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0))
2994 Matrix->Singletons--;
2995 Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] *
2996 Matrix->MarkowitzRow[Col];
2997 if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1))
2998 Matrix->Singletons--;
2999
3000 return pElement;
3001 }
3002
3003
3004
3005
3006
3007
3008
3009
3010 /*
3011 * ZERO PIVOT ENCOUNTERED
3012 *
3013 * This routine is called when a singular matrix is found. It then
3014 * records the current row and column and exits.
3015 *
3016 * >>> Returned:
3017 * The error code spSINGULAR or spZERO_DIAG is returned.
3018 *
3019 * >>> Arguments:
3020 * Matrix <input> (MatrixPtr)
3021 * Pointer to matrix.
3022 * Step <input> (int)
3023 * Index of diagonal that is zero.
3024 */
3025
3026 static int
MatrixIsSingular(Matrix,Step)3027 MatrixIsSingular( Matrix, Step )
3028
3029 MatrixPtr Matrix;
3030 int Step;
3031 {
3032 /* Begin `MatrixIsSingular'. */
3033
3034 Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
3035 Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
3036 return (Matrix->Error = spSINGULAR);
3037 }
3038
3039
3040 static int
ZeroPivot(Matrix,Step)3041 ZeroPivot( Matrix, Step )
3042
3043 MatrixPtr Matrix;
3044 int Step;
3045 {
3046 /* Begin `ZeroPivot'. */
3047
3048 Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
3049 Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
3050 return (Matrix->Error = spZERO_DIAG);
3051 }
3052
3053
3054
3055
3056
3057
3058 #if ANNOTATE == FULL
3059
3060 /*
3061 *
3062 * WRITE STATUS
3063 *
3064 * Write a summary of important variables to standard output.
3065 */
3066
3067 static
WriteStatus(Matrix,Step)3068 WriteStatus( Matrix, Step )
3069
3070 MatrixPtr Matrix;
3071 int Step;
3072 {
3073 int I;
3074
3075 /* Begin `WriteStatus'. */
3076
3077 printf("Step = %1d ", Step);
3078 printf("Pivot found at %1d,%1d using ", Matrix->PivotsOriginalRow,
3079 Matrix->PivotsOriginalCol);
3080 switch(Matrix->PivotSelectionMethod)
3081 { case 's': printf("SearchForSingleton\n"); break;
3082 case 'q': printf("QuicklySearchDiagonal\n"); break;
3083 case 'd': printf("SearchDiagonal\n"); break;
3084 case 'e': printf("SearchEntireMatrix\n"); break;
3085 }
3086
3087 printf("MarkowitzRow = ");
3088 for (I = 1; I <= Matrix->Size; I++)
3089 printf("%2d ", Matrix->MarkowitzRow[I]);
3090 printf("\n");
3091
3092 printf("MarkowitzCol = ");
3093 for (I = 1; I <= Matrix->Size; I++)
3094 printf("%2d ", Matrix->MarkowitzCol[I]);
3095 printf("\n");
3096
3097 printf("MarkowitzProduct = ");
3098 for (I = 1; I <= Matrix->Size; I++)
3099 printf("%2d ", Matrix->MarkowitzProd[I]);
3100 printf("\n");
3101
3102 printf("Singletons = %2d\n", Matrix->Singletons);
3103
3104 printf("IntToExtRowMap = ");
3105 for (I = 1; I <= Matrix->Size; I++)
3106 printf("%2d ", Matrix->IntToExtRowMap[I]);
3107 printf("\n");
3108
3109 printf("IntToExtColMap = ");
3110 for (I = 1; I <= Matrix->Size; I++)
3111 printf("%2d ", Matrix->IntToExtColMap[I]);
3112 printf("\n");
3113
3114 printf("ExtToIntRowMap = ");
3115 for (I = 1; I <= Matrix->ExtSize; I++)
3116 printf("%2d ", Matrix->ExtToIntRowMap[I]);
3117 printf("\n");
3118
3119 printf("ExtToIntColMap = ");
3120 for (I = 1; I <= Matrix->ExtSize; I++)
3121 printf("%2d ", Matrix->ExtToIntColMap[I]);
3122 printf("\n\n");
3123
3124 /* spPrint((char *)Matrix, NO, YES); */
3125
3126 return;
3127
3128 }
3129 #endif /* ANNOTATE == FULL */
3130