2
* SPARSE FORTRAN MODULE
4
* Author: Advising professor:
5
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
8
* This module contains routines that interface Sparse1.3 to a calling
9
* program written in fortran. Almost every externally available Sparse1.3
10
* routine has a counterpart defined in this file, with the name the
11
* same except the `sp' prefix is changed to `sf'. The spADD_ELEMENT
12
* and spADD_QUAD macros are also replaced with the sfAdd1 and sfAdd4
13
* functions defined in this file.
15
* To ease porting this file to different operating systems, the names of
16
* the functions can be easily redefined (search for `Routine Renaming').
17
* A simple example of a FORTRAN program that calls Sparse is included in
18
* this file (search for Example). When interfacing to a FORTRAN program,
19
* the ARRAY_OFFSET option should be set to NO (see spConfig.h).
22
* These interface routines were written by a C programmer who has little
23
* experience with FORTRAN. The routines have had minimal testing.
24
* Any interface between two languages is going to have portability
25
* problems, this one is no exception.
28
* >>> User accessible functions contained in this file:
73
* FORTRAN -- C COMPATIBILITY
75
* Fortran and C data types correspond in the following way:
76
* -- C -- -- FORTRAN --
77
* int INTEGER*4 or INTEGER*2 (machine dependent, usually int*4)
80
* double DOUBLE PRECISION (used by default in preference to float)
82
* The complex number format used by Sparse is compatible with that
83
* used by FORTRAN. C pointers are passed to FORTRAN as longs, they should
84
* not be used in any way in FORTRAN.
90
* Revision and copyright information.
92
* Copyright (c) 1985,86,87,88
93
* by Kenneth S. Kundert and the University of California.
95
* Permission to use, copy, modify, and distribute this software and its
96
* documentation for any purpose and without fee is hereby granted, provided
97
* that the above copyright notice appear in all copies and supporting
98
* documentation and that the authors and the University of California
99
* are properly credited. The authors and the University of California
100
* make no representations as to the suitability of this software for
101
* any purpose. It is provided `as is', without express or implied warranty.
105
static char copyright[] =
106
"Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
107
static char RCSid[] =
108
"@(#)$Header: spFortran.c,v 1.1 88/06/18 11:15:30 kundert Exp $";
117
* >>> Import descriptions:
119
* Macros that customize the sparse matrix routines.
121
* Macros and declarations to be imported by the user.
123
* Matrix type and macro definitions for the sparse matrix routines.
126
#define spINSIDE_SPARSE
127
#include "spConfig.h"
128
#include "spmatrix.h"
139
#include "../machine.h"
140
#define sfCreate C2F(sfcreate)
141
#define sfStripFills C2F(sfstripfills)
142
#define sfDestroy C2F(sfdestroy)
143
#define sfClear C2F(sfclear)
144
#define sfGetElement C2F(sfgetelement)
145
#define sfGetAdmittance C2F(sfgetadmittance)
146
#define sfGetQuad C2F(sfgetquad)
147
#define sfGetOnes C2F(sfgetones)
148
#define sfAdd1Real C2F(sfadd1real)
149
#define sfAdd1Imag C2F(sfadd1imag)
150
#define sfAdd1Complex C2F(sfadd1complex)
151
#define sfAdd4Real C2F(sfadd4real)
152
#define sfAdd4Imag C2F(sfadd4imag)
153
#define sfAdd4Complex C2F(sfadd4complex)
154
#define sfOrderAndFactor C2F(sforderandfactor)
155
#define sfFactor C2F(sffactor)
156
#define sfPartition C2F(sfpartition)
157
#define sfSolve C2F(sfsolve)
158
#define sfSolveTransposed C2F(sfsolvetransposed)
159
#define sfPrint C2F(sfprint)
160
#define sfFileMatrix C2F(sffilematrix)
161
#define sfFileVector C2F(sffilevector)
162
#define sfFileStats C2F(sffilestats)
163
#define sfMNA_Preorder C2F(sfmna_preorder)
164
#define sfScale C2F(sfscale)
165
#define sfMultiply C2F(sfmultiply)
166
#define sfDeterminant C2F(sfdeterminant)
167
#define sfError C2F(sferror)
168
#define sfWhereSingular C2F(sfwheresingular)
169
#define sfGetSize C2F(sfgetsize)
170
#define sfSetReal C2F(sfsetreal)
171
#define sfSetComplex C2F(sfsetcomplex)
172
#define sfFillinCount C2F(sffillincount)
173
#define sfElementCount C2F(sfelementcount)
174
#define sfDeleteRowAndCol C2F(sfdeleterowandcol)
175
#define sfPseudoCondition C2F(sfpseudocondition)
176
#define sfCondition C2F(sfcondition)
177
#define sfNorm C2F(sfnorm)
178
#define sfLargestElement C2F(sflargestelement)
179
#define sfRoundoff C2F(sfroundoff)
184
* Example of a FORTRAN Program Calling Sparse
187
integer matrix, error, sfCreate, sfGetElement, spFactor
189
double precision rhs(4), solution(4)
191
matrix = sfCreate(4,0,error)
192
element(1) = sfGetElement(matrix,1,1)
193
element(2) = sfGetElement(matrix,1,2)
194
element(3) = sfGetElement(matrix,2,1)
195
element(4) = sfGetElement(matrix,2,2)
196
element(5) = sfGetElement(matrix,2,3)
197
element(6) = sfGetElement(matrix,3,2)
198
element(7) = sfGetElement(matrix,3,3)
199
element(8) = sfGetElement(matrix,3,4)
200
element(9) = sfGetElement(matrix,4,3)
201
element(10) = sfGetElement(matrix,4,4)
203
call sfAdd1Real(element(1), 2d0)
204
call sfAdd1Real(element(2), -1d0)
205
call sfAdd1Real(element(3), -1d0)
206
call sfAdd1Real(element(4), 3d0)
207
call sfAdd1Real(element(5), -1d0)
208
call sfAdd1Real(element(6), -1d0)
209
call sfAdd1Real(element(7), 3d0)
210
call sfAdd1Real(element(8), -1d0)
211
call sfAdd1Real(element(9), -1d0)
212
call sfAdd1Real(element(10), 3d0)
213
call sfprint(matrix, .false., .false.)
218
error = sfFactor(matrix)
219
call sfSolve(matrix, rhs, solution)
220
write (6, 10) rhs(1), rhs(2), rhs(3), rhs(4)
235
* Allocates and initializes the data structures associated with a matrix.
237
* >>> Returned: [INTEGER]
238
* A pointer to the matrix is returned cast into an integer. This pointer
239
* is then passed and used by the other matrix routines to refer to a
240
* particular matrix. If an error occurs, the NULL pointer is returned.
243
* Size <input> (long *) [INTEGER]
244
* Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
245
* Complex <input> (int *) [INTEGER or INTEGER*2]
246
* Type of matrix. If Complex is 0 then the matrix is real, otherwise
247
* the matrix will be complex. Note that if the routines are not set up
248
* to handle the type of matrix requested, then a spPANIC error will occur.
249
* Further note that if a matrix will be both real and complex, it must
250
* be specified here as being complex.
251
* Error <output> (int *) [INTEGER or INTEGER*2]
252
* Returns error flag, needed because function spError() will not work
253
* correctly if spCreate() returns NULL.
255
* >>> Possible errors:
258
* Error is cleared in this routine.
262
sfCreate( Size, Complex, Error )
264
int *Size, *Complex, *Error;
266
/* Begin `sfCreate'. */
267
return (long)spCreate(*Size, *Complex, Error );
276
* MATRIX DEALLOCATION
278
* Deallocates pointers and elements of Matrix.
281
* Matrix <input> (long *) [INTEGER]
282
* Pointer to the matrix frame which is to be removed from memory.
290
/* Begin `sfDestroy'. */
291
spDestroy((char *)*Matrix);
303
* STRIP FILL-INS FROM MATRIX
305
* Strips the matrix of all fill-ins.
308
* Matrix <input> (long *) [INTEGER]
309
* Pointer to the matrix to be stripped.
313
sfStripFills( Matrix )
317
/* Begin `sfStripFills'. */
318
spStripFills((char *)*Matrix);
332
* Sets every element of the matrix to zero and clears the error flag.
335
* Matrix <input> (long *) [INTEGER]
336
* Pointer to matrix that is to be cleared.
344
/* Begin `sfClear'. */
345
spClear((char *)*Matrix);
355
* SINGLE ELEMENT ADDITION TO MATRIX BY INDEX
357
* Finds element [Row,Col] and returns a pointer to it. If element is
358
* not found then it is created and spliced into matrix. This routine
359
* is only to be used after spCreate() and before spMNA_Preorder(),
360
* spFactor() or spOrderAndFactor(). Returns a pointer to the
361
* Real portion of a MatrixElement. This pointer is later used by
362
* sfAddxxxxx() to directly access element.
364
* >>> Returns: [INTEGER]
365
* Returns a pointer to the element. This pointer is then used to directly
366
* access the element during successive builds.
369
* Matrix <input> (long *) [INTEGER]
370
* Pointer to the matrix that the element is to be added to.
371
* Row <input> (int *) [INTEGER or INTEGER*2]
372
* Row index for element. Must be in the range of [0..Size] unless
373
* the options EXPANDABLE or TRANSLATE are used. Elements placed in
374
* row zero are discarded. In no case may Row be less than zero.
375
* Col <input> (int *) [INTEGER or INTEGER*2]
376
* Column index for element. Must be in the range of [0..Size] unless
377
* the options EXPANDABLE or TRANSLATE are used. Elements placed in
378
* column zero are discarded. In no case may Col be less than zero.
380
* >>> Possible errors:
382
* Error is not cleared in this routine.
386
sfGetElement( Matrix, Row, Col )
391
/* Begin `sfGetElement'. */
392
return (long)spGetElement((char *)*Matrix, *Row, *Col);
403
* ADDITION OF ADMITTANCE TO MATRIX BY INDEX
405
* Performs same function as sfGetElement except rather than one
406
* element, all four Matrix elements for a floating component are
407
* added. This routine also works if component is grounded. Positive
408
* elements are placed at [Node1,Node2] and [Node2,Node1]. This
409
* routine is only to be used after sfCreate() and before
410
* sfMNA_Preorder(), sfFactor() or sfOrderAndFactor().
412
* >>> Returns: [INTEGER or INTEGER*2]
416
* Matrix <input> (long *) [INTEGER]
417
* Pointer to the matrix that component is to be entered in.
418
* Node1 <input> (int *) [INTEGER or INTEGER*2]
419
* Row and column indices for elements. Must be in the range of [0..Size]
420
* unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
421
* ground node. In no case may Node1 be less than zero.
422
* Node2 <input> (int *) [INTEGER or INTEGER*2]
423
* Row and column indices for elements. Must be in the range of [0..Size]
424
* unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
425
* ground node. In no case may Node2 be less than zero.
426
* Template <output> (long[4]) [INTEGER (4)]
427
* Collection of pointers to four elements that are later used to directly
428
* address elements. User must supply the template, this routine will
433
* Error is not cleared in this routine.
437
sfGetAdmittance( Matrix, Node1, Node2, Template )
439
long *Matrix, Template[4];
442
/* Begin `spGetAdmittance'. */
444
( spGetAdmittance((char *)*Matrix, *Node1, *Node2,
445
(struct spTemplate *)Template )
448
#endif /* QUAD_ELEMENT */
460
* ADDITION OF FOUR ELEMENTS TO MATRIX BY INDEX
462
* Similar to sfGetAdmittance, except that sfGetAdmittance only
463
* handles 2-terminal components, whereas sfGetQuad handles simple
464
* 4-terminals as well. These 4-terminals are simply generalized
465
* 2-terminals with the option of having the sense terminals different
466
* from the source and sink terminals. sfGetQuad adds four
467
* elements to the matrix. Positive elements occur at Row1,Col1
468
* Row2,Col2 while negative elements occur at Row1,Col2 and Row2,Col1.
469
* The routine works fine if any of the rows and columns are zero.
470
* This routine is only to be used after sfCreate() and before
471
* sfMNA_Preorder(), sfFactor() or sfOrderAndFactor()
472
* unless TRANSLATE is set true.
474
* >>> Returns: [INTEGER or INTEGER*2]
478
* Matrix <input> (long *) [INTEGER]
479
* Pointer to the matrix that component is to be entered in.
480
* Row1 <input> (int *) [INTEGER or INTEGER*2]
481
* First row index for elements. Must be in the range of [0..Size]
482
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
483
* ground row. In no case may Row1 be less than zero.
484
* Row2 <input> (int *) [INTEGER or INTEGER*2]
485
* Second row index for elements. Must be in the range of [0..Size]
486
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
487
* ground row. In no case may Row2 be less than zero.
488
* Col1 <input> (int *) [INTEGER or INTEGER*2]
489
* First column index for elements. Must be in the range of [0..Size]
490
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
491
* ground column. In no case may Col1 be less than zero.
492
* Col2 <input> (int *) [INTEGER or INTEGER*2]
493
* Second column index for elements. Must be in the range of [0..Size]
494
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
495
* ground column. In no case may Col2 be less than zero.
496
* Template <output> (long[4]) [INTEGER (4)]
497
* Collection of pointers to four elements that are later used to directly
498
* address elements. User must supply the template, this routine will
503
* Error is not cleared in this routine.
507
sfGetQuad( Matrix, Row1, Row2, Col1, Col2, Template )
509
long *Matrix, Template[4];
510
int *Row1, *Row2, *Col1, *Col2;
512
/* Begin `spGetQuad'. */
514
( spGetQuad( (char *)*Matrix, *Row1, *Row2, *Col1, *Col2,
515
(struct spTemplate *)Template )
518
#endif /* QUAD_ELEMENT */
530
* ADDITION OF FOUR STRUCTURAL ONES TO MATRIX BY INDEX
532
* Performs similar function to sfGetQuad() except this routine is
533
* meant for components that do not have an admittance representation.
535
* The following stamp is used:
541
* >>> Returns: [INTEGER or INTEGER*2]
545
* Matrix <input> (long *) [INTEGER]
546
* Pointer to the matrix that component is to be entered in.
547
* Pos <input> (int *) [INTEGER or INTEGER*2]
548
* See stamp above. Must be in the range of [0..Size]
549
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
550
* ground row. In no case may Pos be less than zero.
551
* Neg <input> (int *) [INTEGER or INTEGER*2]
552
* See stamp above. Must be in the range of [0..Size]
553
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
554
* ground row. In no case may Neg be less than zero.
555
* Eqn <input> (int *) [INTEGER or INTEGER*2]
556
* See stamp above. Must be in the range of [0..Size]
557
* unless the options EXPANDABLE or TRANSLATE are used. Zero is the
558
* ground row. In no case may Eqn be less than zero.
559
* Template <output> (long[4]) [INTEGER (4)]
560
* Collection of pointers to four elements that are later used to directly
561
* address elements. User must supply the template, this routine will
566
* Error is not cleared in this routine.
570
sfGetOnes(Matrix, Pos, Neg, Eqn, Template)
572
long *Matrix, Template[4];
573
int *Pos, *Neg, *Eqn;
575
/* Begin `sfGetOnes'. */
577
( spGetOnes( (char *)*Matrix, *Pos, *Neg, *Eqn,
578
(struct spTemplate *)Template )
581
#endif /* QUAD_ELEMENT */
590
* ADD ELEMENT(S) DIRECTLY TO MATRIX
592
* Adds a value to an element or a set of four element in a matrix.
593
* These elements are referenced by pointer, and so must already have
594
* been created by spGetElement(), spGetAdmittance(), spGetQuad(), or
598
* Element <input> (long *) [INTEGER]
599
* Pointer to the element that is to be added to.
600
* Template <input> (long[4]) [INTEGER (4)]
601
* Pointer to the element that is to be added to.
602
* Real <input> (spREAL *) [REAL or DOUBLE PRECISION]
603
* Real portion of the number to be added to the element.
604
* Imag <input> (spREAL *) [REAL or DOUBLE PRECISION]
605
* Imaginary portion of the number to be added to the element.
609
sfAdd1Real( Element, Real )
614
/* Begin `sfAdd1Real'. */
615
*((RealNumber *)*Element) += *Real;
622
sfAdd1Imag( Element, Imag )
627
/* Begin `sfAdd1Imag'. */
628
*(((RealNumber *)*Element)+1) += *Imag;
633
sfAdd1Complex( Element, Real, Imag )
636
RealNumber *Real, *Imag;
638
/* Begin `sfAdd1Complex'. */
639
*((RealNumber *)*Element) += *Real;
640
*(((RealNumber *)*Element)+1) += *Imag;
642
#endif /* spCOMPLEX */
648
sfAdd4Real( Template, Real )
653
/* Begin `sfAdd4Real'. */
654
*((RealNumber *)Template[0]) += *Real;
655
*((RealNumber *)Template[1]) += *Real;
656
*((RealNumber *)Template[2]) -= *Real;
657
*((RealNumber *)Template[3]) -= *Real;
664
sfAdd4Imag( Template, Imag )
669
/* Begin `sfAdd4Imag'. */
670
*(((RealNumber *)Template[0])+1) += *Imag;
671
*(((RealNumber *)Template[1])+1) += *Imag;
672
*(((RealNumber *)Template[2])+1) -= *Imag;
673
*(((RealNumber *)Template[3])+1) -= *Imag;
678
sfAdd4Complex( Template, Real, Imag )
681
RealNumber *Real, *Imag;
683
/* Begin `sfAdd4Complex'. */
684
*((RealNumber *)Template[0]) += *Real;
685
*((RealNumber *)Template[1]) += *Real;
686
*((RealNumber *)Template[2]) -= *Real;
687
*((RealNumber *)Template[3]) -= *Real;
688
*(((RealNumber *)Template[0])+1) += *Imag;
689
*(((RealNumber *)Template[1])+1) += *Imag;
690
*(((RealNumber *)Template[2])+1) -= *Imag;
691
*(((RealNumber *)Template[3])+1) -= *Imag;
693
#endif /* spCOMPLEX */
694
#endif /* QUAD_ELEMENT */
702
* ORDER AND FACTOR MATRIX
704
* This routine chooses a pivot order for the matrix and factors it
705
* into LU form. It handles both the initial factorization and subsequent
706
* factorizations when a reordering is desired. This is handled in a manner
707
* that is transparent to the user. The routine uses a variation of
708
* Gauss's method where the pivots are associated with L and the
709
* diagonal terms of U are one.
711
* >>> Returned: [INTEGER of INTEGER*2]
712
* The error code is returned. Possible errors are listed below.
715
* Matrix <input> (long *) [INTEGER]
717
* RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
718
* Representative right-hand side vector that is used to determine
719
* pivoting order when the right hand side vector is sparse. If
720
* RHS is a NULL pointer then the RHS vector is assumed to
721
* be full and it is not used when determining the pivoting
723
* RelThreshold <input> (RealNumber *) [REAL or DOUBLE PRECISION]
724
* This number determines what the pivot relative threshold will
725
* be. It should be between zero and one. If it is one then the
726
* pivoting method becomes complete pivoting, which is very slow
727
* and tends to fill up the matrix. If it is set close to zero
728
* the pivoting method becomes strict Markowitz with no
729
* threshold. The pivot threshold is used to eliminate pivot
730
* candidates that would cause excessive element growth if they
731
* were used. Element growth is the cause of roundoff error.
732
* Element growth occurs even in well-conditioned matrices.
733
* Setting the RelThreshold large will reduce element growth and
734
* roundoff error, but setting it too large will cause execution
735
* time to be excessive and will result in a large number of
736
* fill-ins. If this occurs, accuracy can actually be degraded
737
* because of the large number of operations required on the
738
* matrix due to the large number of fill-ins. A good value seems
739
* to be 0.001. The default is chosen by giving a value larger
740
* than one or less than or equal to zero. This value should be
741
* increased and the matrix resolved if growth is found to be
742
* excessive. Changing the pivot threshold does not improve
743
* performance on matrices where growth is low, as is often the
744
* case with ill-conditioned matrices. Once a valid threshold is
745
* given, it becomes the new default. The default value of
746
* RelThreshold was choosen for use with nearly diagonally
747
* dominant matrices such as node- and modified-node admittance
748
* matrices. For these matrices it is usually best to use
749
* diagonal pivoting. For matrices without a strong diagonal, it
750
* is usually best to use a larger threshold, such as 0.01 or
752
* AbsThreshold <input> (RealNumber *) [REAL or DOUBLE PRECISION]
753
* The absolute magnitude an element must have to be considered
754
* as a pivot candidate, except as a last resort. This number
755
* should be set significantly smaller than the smallest diagonal
756
* element that is is expected to be placed in the matrix. If
757
* there is no reasonable prediction for the lower bound on these
758
* elements, then AbsThreshold should be set to zero.
759
* AbsThreshold is used to reduce the possibility of choosing as a
760
* pivot an element that has suffered heavy cancellation and as a
761
* result mainly consists of roundoff error. Once a valid
762
* threshold is given, it becomes the new default.
763
* DiagPivoting <input> (long *) [LOGICAL]
764
* A flag indicating that pivot selection should be confined to the
765
* diagonal if possible. If DiagPivoting is nonzero and if
766
* DIAGONAL_PIVOTING is enabled pivots will be chosen only from
767
* the diagonal unless there are no diagonal elements that satisfy
768
* the threshold criteria. Otherwise, the entire reduced
769
* submatrix is searched when looking for a pivot. The diagonal
770
* pivoting in Sparse is efficient and well refined, while the
771
* off-diagonal pivoting is not. For symmetric and near symmetric
772
* matrices, it is best to use diagonal pivoting because it
773
* results in the best performance when reordering the matrix and
774
* when factoring the matrix without ordering. If there is a
775
* considerable amount of nonsymmetry in the matrix, then
776
* off-diagonal pivoting may result in a better equation ordering
777
* simply because there are more pivot candidates to choose from.
778
* A better ordering results in faster subsequent factorizations.
779
* However, the initial pivot selection process takes considerably
780
* longer for off-diagonal pivoting.
782
* >>> Possible errors:
786
* Error is cleared in this function.
790
sfOrderAndFactor( Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )
792
long *Matrix, *DiagPivoting;
793
RealNumber RHS[], *RelThreshold, *AbsThreshold;
795
/* Begin `sfOrderAndFactor'. */
796
return spOrderAndFactor( (char *)*Matrix, RHS, *RelThreshold,
797
*AbsThreshold, (BOOLEAN)*DiagPivoting );
809
* This routine is the companion routine to spOrderAndFactor().
810
* Unlike sfOrderAndFactor(), sfFactor() cannot change the ordering.
811
* It is also faster than sfOrderAndFactor(). The standard way of
812
* using these two routines is to first use sfOrderAndFactor() for the
813
* initial factorization. For subsequent factorizations, sfFactor()
814
* is used if there is some assurance that little growth will occur
815
* (say for example, that the matrix is diagonally dominant). If
816
* sfFactor() is called for the initial factorization of the matrix,
817
* then sfOrderAndFactor() is automatically called with the default
818
* threshold. This routine uses "row at a time" LU factorization.
819
* Pivots are associated with the lower triangular matrix and the
820
* diagonals of the upper triangular matrix are ones.
822
* >>> Returned: [INTEGER or INTEGER*2]
823
* The error code is returned. Possible errors are listed below.
826
* Matrix <input> (long *) [INTEGER]
829
* >>> Possible errors:
834
* Error is cleared in this function.
842
/* Begin `sfFactor'. */
843
return spFactor((char *)*Matrix);
854
* This routine determines the cost to factor each row using both
855
* direct and indirect addressing and decides, on a row-by-row basis,
856
* which addressing mode is fastest. This information is used in
857
* sfFactor() to speed the factorization.
859
* When factoring a previously ordered matrix using sfFactor(), Sparse
860
* operates on a row-at-a-time basis. For speed, on each step, the
861
* row being updated is copied into a full vector and the operations
862
* are performed on that vector. This can be done one of two ways,
863
* either using direct addressing or indirect addressing. Direct
864
* addressing is fastest when the matrix is relatively dense and
865
* indirect addressing is best when the matrix is quite sparse. The
866
* user selects the type of partition used with Mode. If Mode is set
867
* to spDIRECT_PARTITION, then the all rows are placed in the direct
868
* addressing partition. Similarly, if Mode is set to
869
* spINDIRECT_PARTITION, then the all rows are placed in the indirect
870
* addressing partition. By setting Mode to spAUTO_PARTITION, the
871
* user allows Sparse to select the partition for each row
872
* individually. sfFactor() generally runs faster if Sparse is
873
* allowed to choose its own partitioning, however choosing a
874
* partition is expensive. The time required to choose a partition is
875
* of the same order of the cost to factor the matrix. If you plan to
876
* factor a large number of matrices with the same structure, it is
877
* best to let Sparse choose the partition. Otherwise, you should
878
* choose the partition based on the predicted density of the matrix.
881
* Matrix <input> (long *) [INTEGER]
883
* Mode <input> (int *) [INTEGER or INTEGER*2]
884
* Mode must be one of three special codes: spDIRECT_PARTITION,
885
* spINDIRECT_PARTITION, or spAUTO_PARTITION.
889
sfPartition( Matrix, Mode )
894
/* Begin `sfPartition'. */
895
spPartition((char *)*Matrix, *Mode);
905
* SOLVE MATRIX EQUATION
907
* Performs forward elimination and back substitution to find the
908
* unknown vector from the RHS vector and factored matrix. This
909
* routine assumes that the pivots are associated with the lower
910
* triangular (L) matrix and that the diagonal of the upper triangular
911
* (U) matrix consists of ones. This routine arranges the computation
912
* in different way than is traditionally used in order to exploit the
913
* sparsity of the right-hand side. See the reference in spRevision.
916
* Matrix <input> (long *) [INTEGER]
918
* RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
919
* RHS is the input data array, the right hand side. This data is
920
* undisturbed and may be reused for other solves.
921
* Solution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
922
* Solution is the output data array. This routine is constructed such that
923
* RHS and Solution can be the same array.
924
* iRHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
925
* iRHS is the imaginary portion of the input data array, the right
926
* hand side. This data is undisturbed and may be reused for other solves.
927
* This argument is only necessary if matrix is complex and if
928
* spSEPARATED_COMPLEX_VECTOR is set true.
929
* iSolution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
930
* iSolution is the imaginary portion of the output data array. This
931
* routine is constructed such that iRHS and iSolution can be
932
* the same array. This argument is only necessary if matrix is complex
933
* and if spSEPARATED_COMPLEX_VECTOR is set true.
937
* Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
938
* spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
945
sfSolve( Matrix, RHS, Solution IMAG_VECTORS )
948
RealVector RHS, Solution IMAG_VECTORS;
950
/* Begin `sfSolve'. */
951
spSolve( (char *)*Matrix, RHS, Solution IMAG_VECTORS );
961
* SOLVE TRANSPOSED MATRIX EQUATION
963
* Performs forward elimination and back substitution to find the
964
* unknown vector from the RHS vector and transposed factored
965
* matrix. This routine is useful when performing sensitivity analysis
966
* on a circuit using the adjoint method. This routine assumes that
967
* the pivots are associated with the untransposed lower triangular
968
* (L) matrix and that the diagonal of the untransposed upper
969
* triangular (U) matrix consists of ones.
972
* Matrix <input> (long *) [INTEGER]
974
* RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
975
* RHS is the input data array, the right hand side. This data is
976
* undisturbed and may be reused for other solves.
977
* Solution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
978
* Solution is the output data array. This routine is constructed such that
979
* RHS and Solution can be the same array.
980
* iRHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
981
* iRHS is the imaginary portion of the input data array, the right
982
* hand side. This data is undisturbed and may be reused for other solves.
983
* If spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real, there
984
* is no need to supply this array.
985
* iSolution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)]
986
* iSolution is the imaginary portion of the output data array. This
987
* routine is constructed such that iRHS and iSolution can be
988
* the same array. If spSEPARATED_COMPLEX_VECTOR is set false, or if
989
* matrix is real, there is no need to supply this array.
993
* Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
994
* spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
1001
sfSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS )
1004
RealVector RHS, Solution IMAG_VECTORS;
1006
/* Begin `sfSolveTransposed'. */
1007
spSolveTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS );
1009
#endif /* TRANSPOSE */
1019
* Formats and send the matrix to standard output. Some elementary
1020
* statistics are also output. The matrix is output in a format that is
1021
* readable by people.
1024
* Matrix <input> (long *) [INTEGER]
1025
* Pointer to matrix.
1026
* PrintReordered <input> (long *) [LOGICAL]
1027
* Indicates whether the matrix should be printed out in its original
1028
* form, as input by the user, or whether it should be printed in its
1029
* reordered form, as used by the matrix routines. A zero indicates that
1030
* the matrix should be printed as inputed, a one indicates that it
1031
* should be printed reordered.
1032
* Data <input> (long *) [LOGICAL]
1033
* Boolean flag that when false indicates that output should be
1034
* compressed such that only the existence of an element should be
1035
* indicated rather than giving the actual value. Thus 10 times as many
1036
* can be printed on a row. A zero signifies that the matrix should
1037
* be printed compressed. A one indicates that the matrix should be
1038
* printed in all its glory.
1039
* Header <input> (long *) [LOGICAL]
1040
* Flag indicating that extra information such as the row and column
1041
* numbers should be printed.
1045
sfPrint( Matrix, Data, PrintReordered, Header )
1047
long *Matrix, *PrintReordered, *Data, *Header;
1049
/* Begin `sfPrint'. */
1050
spPrint( (char *)*Matrix, (int)*PrintReordered, (int)*Data, (int)*Header );
1052
#endif /* DOCUMENTATION */
1061
* OUTPUT MATRIX TO FILE
1063
* Writes matrix to file in format suitable to be read back in by the
1064
* matrix test program. Data is sent to a file with a fixed name because
1065
* it is impossible to pass strings from FORTRAN to C in a manner that is
1069
* One is returned if routine was successful, otherwise zero is returned.
1070
* The calling function can query errno (the system global error variable)
1071
* as to the reason why this routine failed.
1073
* >>> Arguments: [LOGICAL]
1074
* Matrix <input> (long *) [INTEGER]
1075
* Pointer to matrix.
1076
* Reordered <input> (long *) [LOGICAL]
1077
* Specifies whether matrix should be output in reordered form,
1078
* or in original order.
1079
* Data <input> (long *) [LOGICAL]
1080
* Indicates that the element values should be output along with
1081
* the indices for each element. This parameter must be true if
1082
* matrix is to be read by the sparse test program.
1083
* Header <input> (long *) [LOGICAL]
1084
* Indicates that header is desired. This parameter must be true if
1085
* matrix is to be read by the sparse test program.
1087
#define MATRIX_FILE_NAME "spMatrix"
1088
#define STATS_FILE_NAME "spStats"
1091
sfFileMatrix( Matrix, Reordered, Data, Header )
1093
long *Matrix, *Reordered, *Data, *Header;
1095
/* Begin `sfFileMatrix'. */
1096
return spFileMatrix( (char *)*Matrix, MATRIX_FILE_NAME, "",
1097
(int)*Reordered, (int)*Data, (int)*Header );
1099
#endif /* DOCUMENTATION */
1105
* OUTPUT SOURCE VECTOR TO FILE
1107
* Writes vector to file in format suitable to be read back in by the
1108
* matrix test program. This routine should be executed after the function
1112
* One is returned if routine was successful, otherwise zero is returned.
1113
* The calling function can query errno (the system global error variable)
1114
* as to the reason why this routine failed.
1117
* Matrix <input> (long *)
1118
* Pointer to matrix.
1119
* RHS <input> (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)]
1120
* Right-hand side vector. This is only the real portion if
1121
* spSEPARATED_COMPLEX_VECTORS is true.
1122
* iRHS <input> (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)]
1123
* Right-hand side vector, imaginary portion. Not necessary if matrix
1124
* is real or if spSEPARATED_COMPLEX_VECTORS is set false.
1128
sfFileVector( Matrix, RHS IMAG_RHS )
1131
RealVector RHS IMAG_RHS;
1133
/* Begin `sfFileVector'. */
1134
return spFileVector( (char *)*Matrix, MATRIX_FILE_NAME, RHS IMAG_RHS );
1136
#endif /* DOCUMENTATION */
1146
* OUTPUT STATISTICS TO FILE
1148
* Writes useful information concerning the matrix to a file. Should be
1149
* executed after the matrix is factored.
1151
* >>> Returns: [LOGICAL]
1152
* One is returned if routine was successful, otherwise zero is returned.
1153
* The calling function can query errno (the system global error variable)
1154
* as to the reason why this routine failed.
1157
* Matrix <input> (long *) [INTEGER]
1158
* Pointer to matrix.
1162
sfFileStats( Matrix )
1166
/* Begin `sfFileStats'. */
1167
return spFileStats( (char *)*Matrix, STATS_FILE_NAME, "" );
1169
#endif /* DOCUMENTATION */
1176
* PREORDER MODIFIED NODE ADMITTANCE MATRIX TO REMOVE ZEROS FROM DIAGONAL
1178
* This routine massages modified node admittance matrices to remove
1179
* zeros from the diagonal. It takes advantage of the fact that the
1180
* row and column associated with a zero diagonal usually have
1181
* structural ones placed symmetricly. This routine should be used
1182
* only on modified node admittance matrices and should be executed
1183
* after the matrix has been built but before the factorization
1184
* begins. It should be executed for the initial factorization only
1185
* and should be executed before the rows have been linked. Thus it
1186
* should be run before using spScale(), spMultiply(),
1187
* spDeleteRowAndCol(), or spNorm().
1189
* This routine exploits the fact that the structural one are placed
1190
* in the matrix in symmetric twins. For example, the stamps for
1191
* grounded and a floating voltage sources are
1192
* grounded: floating:
1193
* [ x x 1 ] [ x x 1 ]
1194
* [ x x ] [ x x -1 ]
1196
* Notice for the grounded source, there is one set of twins, and for
1197
* the grounded, there are two sets. We remove the zero from the diagonal
1198
* by swapping the rows associated with a set of twins. For example:
1199
* grounded: floating 1: floating 2:
1200
* [ 1 ] [ 1 -1 ] [ x x 1 ]
1201
* [ x x ] [ x x -1 ] [ 1 -1 ]
1202
* [ x x 1 ] [ x x 1 ] [ x x -1 ]
1204
* It is important to deal with any zero diagonals that only have one
1205
* set of twins before dealing with those that have more than one because
1206
* swapping row destroys the symmetry of any twins in the rows being
1207
* swapped, which may limit future moves. Consider
1212
* There is one set of twins for diagonal 4 and two for diagonal3.
1213
* Dealing with diagonal for first requires swapping rows 2 and 4.
1218
* We can now deal with diagonal 3 by swapping rows 1 and 3.
1223
* And we are done, there are no zeros left on the diagonal. However, if
1224
* we originally dealt with diagonal 3 first, we could swap rows 2 and 3
1229
* Diagonal 4 no longer has a symmetric twin and we cannot continue.
1231
* So we always take care of lone twins first. When none remain, we
1232
* choose arbitrarily a set of twins for a diagonal with more than one set
1233
* and swap the rows corresponding to that twin. We then deal with any
1234
* lone twins that were created and repeat the procedure until no
1235
* zero diagonals with symmetric twins remain.
1237
* In this particular implementation, columns are swapped rather than rows.
1238
* The algorithm used in this function was developed by Ken Kundert and
1242
* Matrix <input> (long *) [INTEGER]
1243
* Pointer to the matrix to be preordered.
1247
sfMNA_Preorder( Matrix )
1251
/* Begin `sfMNA_Preorder'. */
1252
spMNA_Preorder( (char *)*Matrix );
1254
#endif /* MODIFIED_NODAL */
1265
* This function scales the matrix to enhance the possibility of
1266
* finding a good pivoting order. Note that scaling enhances accuracy
1267
* of the solution only if it affects the pivoting order, so it makes
1268
* no sense to scale the matrix before spFactor(). If scaling is
1269
* desired it should be done before spOrderAndFactor(). There
1270
* are several things to take into account when choosing the scale
1271
* factors. First, the scale factors are directly multiplied against
1272
* the elements in the matrix. To prevent roundoff, each scale factor
1273
* should be equal to an integer power of the number base of the
1274
* machine. Since most machines operate in base two, scale factors
1275
* should be a power of two. Second, the matrix should be scaled such
1276
* that the matrix of element uncertainties is equilibrated. Third,
1277
* this function multiplies the scale factors by the elements, so if
1278
* one row tends to have uncertainties 1000 times smaller than the
1279
* other rows, then its scale factor should be 1024, not 1/1024.
1280
* Fourth, to save time, this function does not scale rows or columns
1281
* if their scale factors are equal to one. Thus, the scale factors
1282
* should be normalized to the most common scale factor. Rows and
1283
* columns should be normalized separately. For example, if the size
1284
* of the matrix is 100 and 10 rows tend to have uncertainties near
1285
* 1e-6 and the remaining 90 have uncertainties near 1e-12, then the
1286
* scale factor for the 10 should be 1/1,048,576 and the scale factors
1287
* for the remaining 90 should be 1. Fifth, since this routine
1288
* directly operates on the matrix, it is necessary to apply the scale
1289
* factors to the RHS and Solution vectors. It may be easier to
1290
* simply use spOrderAndFactor() on a scaled matrix to choose the
1291
* pivoting order, and then throw away the matrix. Subsequent
1292
* factorizations, performed with spFactor(), will not need to have
1293
* the RHS and Solution vectors descaled. Lastly, this function
1294
* should not be executed before the function spMNA_Preorder.
1297
* Matrix <input> (long *) [INTEGER]
1298
* Pointer to the matrix to be scaled.
1299
* SolutionScaleFactors <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1300
* The array of Solution scale factors. These factors scale the columns.
1301
* All scale factors are real valued.
1302
* RHS_ScaleFactors <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1303
* The array of RHS scale factors. These factors scale the rows.
1304
* All scale factors are real valued.
1308
sfScale( Matrix, RHS_ScaleFactors, SolutionScaleFactors )
1311
RealVector RHS_ScaleFactors, SolutionScaleFactors;
1313
/* Begin `sfScale'. */
1314
spScale( (char *)*Matrix, RHS_ScaleFactors, SolutionScaleFactors );
1316
#endif /* SCALING */
1325
* MATRIX MULTIPLICATION
1327
* Multiplies matrix by solution vector to find source vector.
1328
* Assumes matrix has not been factored. This routine can be used
1329
* as a test to see if solutions are correct. It should not be used
1330
* before PreorderFoModifiedNodal().
1333
* Matrix <input> (long *) [INTEGER]
1334
* Pointer to the matrix.
1335
* RHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1336
* RHS is the right hand side. This is what is being solved for.
1337
* Solution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1338
* Solution is the vector being multiplied by the matrix.
1339
* iRHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1340
* iRHS is the imaginary portion of the right hand side. This is
1341
* what is being solved for. This is only necessary if the matrix is
1342
* complex and spSEPARATED_COMPLEX_VECTORS is true.
1343
* iSolution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1344
* iSolution is the imaginary portion of the vector being multiplied
1345
* by the matrix. This is only necessary if the matrix is
1346
* complex and spSEPARATED_COMPLEX_VECTORS is true.
1348
* >>> Obscure Macros
1350
* Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
1351
* spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
1356
sfMultiply( Matrix, RHS, Solution IMAG_VECTORS )
1359
RealVector Solution, RHS IMAG_VECTORS;
1361
/* Begin `sfMultiply'. */
1362
spMultiply( (char *)*Matrix, RHS, Solution IMAG_VECTORS );
1364
#endif /* MULTIPLICATION */
1371
#if MULTIPLICATION AND TRANSPOSE
1373
* TRANSPOSED MATRIX MULTIPLICATION
1375
* Multiplies transposed matrix by solution vector to find source vector.
1376
* Assumes matrix has not been factored. This routine can be used
1377
* as a test to see if solutions are correct. It should not be used
1378
* before PreorderFoModifiedNodal().
1381
* Matrix <input> (long *) [INTEGER]
1382
* Pointer to the matrix.
1383
* RHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1384
* RHS is the right hand side. This is what is being solved for.
1385
* Solution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1386
* Solution is the vector being multiplied by the matrix.
1387
* iRHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1388
* iRHS is the imaginary portion of the right hand side. This is
1389
* what is being solved for. This is only necessary if the matrix is
1390
* complex and spSEPARATED_COMPLEX_VECTORS is true.
1391
* iSolution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)]
1392
* iSolution is the imaginary portion of the vector being multiplied
1393
* by the matrix. This is only necessary if the matrix is
1394
* complex and spSEPARATED_COMPLEX_VECTORS is true.
1396
* >>> Obscure Macros
1398
* Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
1399
* spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
1404
sfMultTransposed( Matrix, RHS, Solution IMAG_VECTORS )
1407
RealVector Solution, RHS IMAG_VECTORS;
1409
/* Begin `sfMultTransposed'. */
1410
spMultTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS );
1412
#endif /* MULTIPLICATION AND TRANSPOSE */
1422
* CALCULATE DETERMINANT
1424
* This routine in capable of calculating the determinant of the
1425
* matrix once the LU factorization has been performed. Hence, only
1426
* use this routine after spFactor() and before spClear().
1427
* The determinant equals the product of all the diagonal elements of
1428
* the lower triangular matrix L, except that this product may need
1429
* negating. Whether the product or the negative product equals the
1430
* determinant is determined by the number of row and column
1431
* interchanges performed. Note that the determinants of matrices can
1432
* be very large or very small. On large matrices, the determinant
1433
* can be far larger or smaller than can be represented by a floating
1434
* point number. For this reason the determinant is scaled to a
1435
* reasonable value and the logarithm of the scale factor is returned.
1438
* Matrix <input> (long *) [INTEGER]
1439
* A pointer to the matrix for which the determinant is desired.
1440
* pExponent <output> (int *) [INTEGER or INTEGER*2]
1441
* The logarithm base 10 of the scale factor for the determinant. To
1443
* the actual determinant, Exponent should be added to the exponent of
1445
* pDeterminant <output> (RealNumber *) [REAL or DOUBLE PRECISION]
1446
* The real portion of the determinant. This number is scaled to be
1447
* greater than or equal to 1.0 and less than 10.0.
1448
* piDeterminant <output> (RealNumber *) [REAL or DOUBLE PRECISION]
1449
* The imaginary portion of the determinant. When the matrix is real
1450
* this pointer need not be supplied, nothing will be returned. This
1451
* number is scaled to be greater than or equal to 1.0 and less than 10.0.
1457
sfDeterminant( Matrix, pExponent, pDeterminant, piDeterminant )
1460
RealNumber *pDeterminant, *piDeterminant;
1463
/* Begin `sfDeterminant'. */
1464
spDeterminant( (char *)*Matrix, pExponent, pDeterminant, piDeterminant );
1467
#else /* spCOMPLEX */
1470
sfDeterminant( Matrix, pExponent, pDeterminant )
1473
RealNumber *pDeterminant;
1476
/* Begin `sfDeterminant'. */
1477
spDeterminant( (char *)*Matrix, pExponent, pDeterminant );
1479
#endif /* spCOMPLEX */
1480
#endif /* DETERMINANT */
1488
* RETURN MATRIX ERROR STATUS
1490
* This function is used to determine the error status of the given matrix.
1492
* >>> Returned: [INTEGER or INTEGER*2]
1493
* The error status of the given matrix.
1496
* Matrix <input> (long *) [INTEGER]
1497
* The matrix for which the error status is desired.
1505
/* Begin `sfError'. */
1506
return spError( (char *)*Matrix );
1515
* WHERE IS MATRIX SINGULAR
1517
* This function returns the row and column number where the matrix was
1518
* detected as singular or where a zero was detected on the diagonal.
1521
* Matrix <input> (long *) [INTEGER]
1522
* The matrix for which the error status is desired.
1523
* pRow <output> (int *) [INTEGER or INTEGER*2]
1525
* pCol <output> (int *) [INTEGER or INTEGER*2]
1526
* The column number.
1530
sfWhereSingular( Matrix, Row, Col )
1535
/* Begin `sfWhereSingular'. */
1536
spWhereSingular( (char *)*Matrix, Row, Col );
1546
* Returns the size of the matrix. Either the internal or external size of
1547
* the matrix is returned.
1550
* Matrix <input> (long *) [INTEGER]
1551
* Pointer to matrix.
1552
* External <input> (BOOLEAN) [LOGICAL]
1553
* If External is set true, the external size , i.e., the value of the
1554
* largest external row or column number encountered is returned.
1555
* Otherwise the true size of the matrix is returned. These two sizes
1556
* may differ if the TRANSLATE option is set true.
1560
sfGetSize( Matrix, External )
1562
long *Matrix, *External;
1564
/* Begin `sfGetSize'. */
1565
return spGetSize( (char *)*Matrix, (BOOLEAN)*External );
1576
* SET MATRIX COMPLEX OR REAL
1578
* Forces matrix to be either real or complex.
1581
* Matrix <input> (long *) [INTEGER]
1582
* Pointer to matrix.
1590
/* Begin `sfSetReal'. */
1591
spSetReal( (char *)*Matrix );
1596
sfSetComplex( Matrix )
1600
/* Begin `sfSetComplex'. */
1601
spSetComplex( (char *)*Matrix );
1613
* ELEMENT OR FILL-IN COUNT
1615
* Two functions used to return simple statistics. Either the number
1616
* of total elements, or the number of fill-ins can be returned.
1619
* Matrix <input> (long *) [INTEGER]
1620
* Pointer to matrix.
1624
sfFillinCount( Matrix )
1628
/* Begin `sfFillinCount'. */
1629
return spFillinCount( (char *)*Matrix );
1634
sfElementCount( Matrix )
1638
/* Begin `sfElementCount'. */
1639
return spElementCount( (char *)*Matrix );
1647
#if TRANSLATE AND DELETE
1650
* DELETE A ROW AND COLUMN FROM THE MATRIX
1652
* Deletes a row and a column from a matrix.
1654
* Sparse will abort if an attempt is made to delete a row or column that
1658
* Matrix <input> (long *) [INTEGER]
1659
* Pointer to the matrix in which the row and column are to be deleted.
1660
* Row <input> (int) [INTEGER or INTEGER*2]
1661
* Row to be deleted.
1662
* Col <input> (int) [INTEGER or INTEGER*2]
1663
* Column to be deleted.
1667
sfDeleteRowAndCol( Matrix, Row, Col )
1672
/* Begin `sfDeleteRowAndCol'. */
1673
spDeleteRowAndCol( (char *)*Matrix, *Row, *Col );
1684
* CALCULATE PSEUDOCONDITION
1686
* Computes the magnitude of the ratio of the largest to the smallest
1687
* pivots. This quantity is an indicator of ill-conditioning in the
1688
* matrix. If this ratio is large, and if the matrix is scaled such
1689
* that uncertainties in the RHS and the matrix entries are
1690
* equilibrated, then the matrix is ill-conditioned. However, a small
1691
* ratio does not necessarily imply that the matrix is
1692
* well-conditioned. This routine must only be used after a matrix
1693
* has been factored by sfOrderAndFactor() or sfFactor() and before it
1694
* is cleared by sfClear() or spInitialize(). The pseudocondition is faster
1695
* to compute than the condition number calculated by sfCondition(), but
1696
* is not as informative.
1698
* >>> Returns: [REAL or DOUBLE PRECISION]
1699
* The magnitude of the ratio of the largest to smallest pivot used during
1700
* previous factorization. If the matrix was singular, zero is returned.
1703
* Matrix <input> (long *)
1704
* Pointer to the matrix.
1708
sfPseudoCondition( Matrix )
1712
/* Begin `sfPseudoCondition'. */
1713
return spPseudoCondition( (char *)Matrix );
1726
* ESTIMATE CONDITION NUMBER
1728
* Computes an estimate of the condition number using a variation on
1729
* the LINPACK condition number estimation algorithm. This quantity is
1730
* an indicator of ill-conditioning in the matrix. To avoid problems
1731
* with overflow, the reciprocal of the condition number is returned.
1732
* If this number is small, and if the matrix is scaled such that
1733
* uncertainties in the RHS and the matrix entries are equilibrated,
1734
* then the matrix is ill-conditioned. If the this number is near
1735
* one, the matrix is well conditioned. This routine must only be
1736
* used after a matrix has been factored by sfOrderAndFactor() or
1737
* sfFactor() and before it is cleared by sfClear() or spInitialize().
1739
* Unlike the LINPACK condition number estimator, this routines
1740
* returns the L infinity condition number. This is an artifact of
1741
* Sparse placing ones on the diagonal of the upper triangular matrix
1742
* rather than the lower. This difference should be of no importance.
1745
* A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate
1746
* for the condition number of a matrix. SIAM Journal on Numerical
1747
* Analysis. Vol. 16, No. 2, pages 368-375, April 1979.
1749
* J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK
1750
* User's Guide. SIAM, 1979.
1752
* Roger G. Grimes, John G. Lewis. Condition number estimation for
1753
* sparse matrices. SIAM Journal on Scientific and Statistical
1754
* Computing. Vol. 2, No. 4, pages 384-388, December 1981.
1756
* Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM
1757
* Journal on Scientific and Statistical Computing. Vol. 1, No. 2,
1758
* pages 205-209, June 1980.
1760
* >>> Returns: [REAL or DOUBLE PRECISION]
1761
* The reciprocal of the condition number. If the matrix was singular,
1765
* eMatrix <input> (long *)
1766
* Pointer to the matrix.
1767
* NormOfMatrix <input> (RealNumber *) [REAL or DOUBLE PRECISION]
1768
* The L-infinity norm of the unfactored matrix as computed by
1770
* pError <output> (int *) [INTEGER or INTEGER*2]
1771
* Used to return error code.
1773
* >>> Possible errors:
1779
sfCondition( Matrix, NormOfMatrix, pError )
1782
RealNumber *NormOfMatrix;
1785
/* Begin `sfCondition'. */
1786
return spCondition( (char *)*Matrix, *NormOfMatrix, pError );
1794
* L-INFINITY MATRIX NORM
1796
* Computes the L-infinity norm of an unfactored matrix. It is a fatal
1797
* error to pass this routine a factored matrix.
1799
* One difficulty is that the rows may not be linked.
1801
* >>> Returns: [REAL or DOUBLE PRECISION]
1802
* The largest absolute row sum of matrix.
1805
* Matrix <input> (long *)
1806
* Pointer to the matrix.
1814
/* Begin `sfNorm'. */
1815
return spNorm( (char *)*Matrix );
1817
#endif /* CONDITION */
1826
* STABILITY OF FACTORIZATION
1828
* The following routines are used to gauge the stability of a
1829
* factorization. If the factorization is determined to be too unstable,
1830
* then the matrix should be reordered. The routines compute quantities
1831
* that are needed in the computation of a bound on the error attributed
1832
* to any one element in the matrix during the factorization. In other
1833
* words, there is a matrix E = [e_ij] of error terms such that A+E = LU.
1834
* This routine finds a bound on |e_ij|. Erisman & Reid [1] showed that
1835
* |e_ij| < 3.01 u rho m_ij, where u is the machine rounding unit,
1836
* rho = max a_ij where the max is taken over every row i, column j, and
1837
* step k, and m_ij is the number of multiplications required in the
1838
* computation of l_ij if i > j or u_ij otherwise. Barlow [2] showed that
1839
* rho < max_i || l_i ||_p max_j || u_j ||_q where 1/p + 1/q = 1.
1841
* The first routine finds the magnitude on the largest element in the
1842
* matrix. If the matrix has not yet been factored, the largest
1843
* element is found by direct search. If the matrix is factored, a
1844
* bound on the largest element in any of the reduced submatrices is
1845
* computed using Barlow with p = oo and q = 1. The ratio of these
1846
* two numbers is the growth, which can be used to determine if the
1847
* pivoting order is adequate. A large growth implies that
1848
* considerable error has been made in the factorization and that it
1849
* is probably a good idea to reorder the matrix. If a large growth
1850
* in encountered after using spFactor(), reconstruct the matrix and
1851
* refactor using spOrderAndFactor(). If a large growth is
1852
* encountered after using spOrderAndFactor(), refactor using
1853
* spOrderAndFactor() with the pivot threshold increased, say to 0.1.
1855
* Using only the size of the matrix as an upper bound on m_ij and
1856
* Barlow's bound, the user can estimate the size of the matrix error
1857
* terms e_ij using the bound of Erisman and Reid. The second routine
1858
* computes a tighter bound (with more work) based on work by Gear
1859
* [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the
1860
* threshold and c is the maximum number of off-diagonal elements in
1861
* any row of L. The expensive part of computing this bound is
1862
* determining the maximum number of off-diagonals in L, which changes
1863
* only when the order of the matrix changes. This number is computed
1864
* and saved, and only recomputed if the matrix is reordered.
1866
* [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the
1867
* triangular factorization of a sparse matrix. Numerische
1868
* Mathematik. Vol. 22, No. 3, 1974, pp 183-186.
1870
* [2] J. L. Barlow. A note on monitoring the stability of triangular
1871
* decomposition of sparse matrices. "SIAM Journal of Scientific
1872
* and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168.
1874
* [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse
1875
* Matrices." Oxford 1986. pp 99.
1879
* LARGEST ELEMENT IN MATRIX
1881
* >>> Returns: [REAL or DOUBLE PRECISION]
1882
* If matrix is not factored, returns the magnitude of the largest element in
1883
* the matrix. If the matrix is factored, a bound on the magnitude of the
1884
* largest element in any of the reduced submatrices is returned.
1887
* Matrix <input> (long *) [INTEGER]
1888
* Pointer to the matrix.
1892
sfLargestElement( Matrix )
1896
/* Begin `sfLargestElement'. */
1897
return spLargestElement( (char *)Matrix );
1904
* MATRIX ROUNDOFF ERROR
1906
* >>> Returns: [REAL or DOUBLE PRECISION]
1907
* Returns a bound on the magnitude of the largest element in E = A - LU.
1910
* Matrix <input> (long *) [INTEGER]
1911
* Pointer to the matrix.
1912
* Rho <input> (RealNumber *) [REAL or DOUBLE PRECISION]
1913
* The bound on the magnitude of the largest element in any of the
1914
* reduced submatrices. This is the number computed by the function
1915
* spLargestElement() when given a factored matrix. If this number is
1916
* negative, the bound will be computed automatically.
1920
sfRoundoff( Matrix, Rho )
1925
/* Begin `sfRoundoff'. */
1926
return spRoundoff( (char *)*Matrix, *Rho );
1930
#endif /* FORTRAN */
added
removed
routines/sparse/spFortran.c
