1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 Sparse User's Guide 20 21 A Sparse Linear Equation Solver 22 23 24 Version 1.3a 25 26 1 April 1988 27 28 29 30 31 32 Kenneth S. Kundert 33 Alberto Sangiovanni-Vincentelli 34 35 36 37 38 39 40 Department of 41 Electrical Engineering and Computer Sciences 42 University of California, Berkeley 43 Berkeley, Calif. 94720 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 June 23, 1988 67 68 69 70 71 72 73 74 751: INTRODUCTION 76 77 Sparse1.3 is a flexible package of subroutines written in C used to 78quickly and accurately solve large sparse systems of linear equations. The 79package is able to handle arbitrary real and complex square matrix equa- 80tions. Besides being able to solve linear systems, it is also able to 81quickly solve transposed systems, find determinants, and estimate errors 82due to ill-conditioning in the system of equations and instability in the 83computations. Sparse also provides a test program that is able read matrix 84equations from a file, solve them, and print useful information about the 85equation and its solution. 86 87 Sparse1.3 is generally as fast or faster than other popular sparse 88matrix packages when solving many matrices of similar structure. Sparse 89does not require or assume symmetry and is able to perform numerical pivot- 90ing to avoid unnecessary error in the solution. It handles its own memory 91allocation, which allows the user to forgo the hassle of providing adequate 92memory. It also has a natural, flexible, and efficient interface to the 93calling program. 94 95 Sparse was originally written for use in circuit simulators and is 96particularly apt at handling node- and modified-node admittance matrices. 97The systems of linear generated in a circuit simulator stem from solving 98large systems of nonlinear equations using Newton's method and integrating 99large stiff systems of ordinary differential equations. However, Sparse is 100also suitable for other uses, one in particular is solving the very large 101systems of linear equations resulting from the numerical solution of par- 102tial differential equations. 103 104 1051.1: Features of Sparse1.3 106 107 Beyond the basic capability of being able to create, factor and solve 108systems of equations, this package features several other capabilities that 109enhance its utility. These features are: 110 111o Ability to handle both real and complex systems of equations. Both 112 types may resident and active at the same time. In fact, the same 113 matrix may alternate between being real and complex. 114 115o Ability to quickly solve the transposed system. This feature is use- 116 ful when computing the sensitivity of a circuit using the adjoint 117 method. 118 119o Memory for elements in the matrix is allocated dynamically, so the 120 size of the matrix is only limited by the amount of memory available 121 to Sparse and the range of the integer data type, which is used to 122 hold matrix indices. 123 124o Ability to efficiently compute the condition number of the matrix and 125 an a posteriori estimate of the error caused by growth in the size of 126 the elements during the factorization. 127 128o Much of the matrix initialization can be performed by Sparse, 129 130 131 132 June 23, 1988 133 134 135 136 137 138 - 2 - 139 140 141 providing advantages in speed and simplified coding of the calling 142 program. 143 144o Ability to preorder modified node admittance matrices to enhance accu- 145 racy and speed. 146 147o Ability to exploit sparsity in the right-hand side vector to reduce 148 unnecessary computation. 149 150o Ability to scale matrices prior to factoring to reduce uncertainty in 151 the solution. 152 153o The ability to create and build a matrix without knowing its final 154 size. 155 156o The ability to add elements, and rows and columns, to a matrix after 157 the matrix has been reordered. 158 159o The ability to delete rows and columns from a matrix. 160 161o The ability to strip the fill-ins from a matrix. This can improve the 162 efficiency of a subsequent reordering. 163 164o The ability to handle matrices that have rows and columns missing from 165 their input description. 166 167o Ability to output the matrix in forms readable by either by people or 168 by the Sparse package. Basic statistics on the matrix can also be 169 output. 170 171o By default all arithmetic operations and number storage use double 172 precision. Thus, Sparse usually gives accurate results, even on 173 highly ill-conditioned systems. If so desired, Sparse can be easily 174 configured to use single precision arithmetic. 175 176 1771.2: Enhancements of Sparse1.3 over Sparse1.2 178 179 Most notable of the enhancements provided by Sparse1.3 is that it is 180considerably faster on dense matrices. Also, external names have been made 181unique to 7 characters and the Sparse prefix sp has been prepended to all 182externally accessible names to avoid conflicts. In addition, a routine 183that efficiently estimates the condition number of a matrix has been added 184and the code that estimates the growth in the factorization has been split 185off from the actual factorization so that it is computed only when needed. 186 187 It is now possible for the user program to store information in the 188matrix elements. It is also possible to provide a subroutine to Sparse 189that uses that information to initialize the matrix. This can greatly sim- 190plify the user's code. 191 192 Though the interface between Sparse1.3 and the calling program has 193changed considerable from previous version of Sparse, it is possible to 194compile additional code that provides backward compatibility to Sparse1.2 195 196 197 198 June 23, 1988 199 200 201 202 203 204 - 3 - 205 206 207at the expense of a slight loss of efficiency by setting a compiler option. 208Sparse1.3 now also has an FORTRAN interface. Routines written in FORTRAN 209can access almost all of the features Sparse1.3. 210 211 2121.3: Copyright Information 213 214 Sparse1.3 has been copyrighted. Permission to use, copy, modify, and 215distribute this software and its documentation for any purpose and without 216fee is hereby granted, provided that the copyright notice appear in all 217copies, and Sparse and the University of California, Berkeley are refer- 218enced in all documentation for the program or product in which Sparse is to 219be installed. The authors and the University of California make no 220representations as to the suitability of the software for any purpose. It 221is provided `as is', without express or implied warranty. 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 June 23, 1988 265 266 267 268 269 270 - 4 - 271 272 2732: PRIMER 274 2752.1: Solving Matrix Equations 276 277 Sparse contains a collection of C subprograms that can be used to 278solve linear algebraic systems of equations. These systems are of the 279form: 280 281 Ax = b 282where A is an nxn matrix, x is the vector of n unknowns and b is the vector 283of n right-hand side terms. Through out this package A is denoted Matrix, 284x is denoted Solution and b is denoted RHS (for right-hand side). The sys- 285tem is solved using LU factorization, so the actual solution process is 286broken into two steps, the factorization or decomposition of the matrix, 287performed by spFactor(), and the forward and backward substitution, per- 288formed by spSolve(). spFactor() factors the given matrix into upper and 289lower triangular matrices independent of the right-hand side. Once this is 290done, the solution vector can be determined efficiently for any number of 291right-hand sides without refactoring the matrix. 292 293 This package exploits the fact that large matrices usually are sparse 294by not storing or operating on elements in the matrix that are zero. Stor- 295ing zero elements is avoided by organizing the matrix into an orthogonal 296linked-list. Thus, to access an element if only its indices are known 297requires stepping through the list, which is slow. This function is per- 298formed by the routine spGetElement(). It is used to initially enter data 299into a matrix and to build the linked-list. Because it is common to 300repeatedly solve matrices with identical zero/nonzero structure, it is pos- 301sible to reuse the linked-list. Thus, the linked list is left in memory 302and the element values are simply cleared by spClear() before the linked- 303list is reused. To speed the entering of the element values into succes- 304sive matrices, spGetElement() returns a pointer to the element in the 305matrix. This pointer can then be used to place data directly into the 306matrix without having to traverse through the linked-list. 307 308 The order in which the rows and columns of the matrix are factored is 309very important. It directly affects the amount of time required for the 310factorization and the forward and backward substitution. It also affects 311the accuracy of the result. The process of choosing this order is time 312consuming, but fortunately it usually only has to be done once for each 313particular matrix structure encountered. When a matrix with a new 314zero/nonzero structure is to be factored, it is done by using spOr- 315derAndFactor(). Subsequent matrices of the same structure are factored 316with spFactor(). The latter routine does not have the ability to reorder 317matrix, but it is considerably faster. It may be that a order chosen may 318be unsuitable for subsequent factorizations. If this is known to be true a 319priori, it is possible to use spOrderAndFactor() for the subsequent factor- 320izations, with a noticeable speed penalty. spOrderAndFactor() monitors the 321numerical stability of the factorization and will modify an existing order- 322ing to maintain stability. Otherwise, an a posteriori measure of the 323numerical stability of the factorization can be computed, and the matrix 324reordered if necessary. 325 326 The Sparse routines allow several matrices of different structures to 327 328 329 330 June 23, 1988 331 332 333 334 335 336 - 5 - 337 338 339be resident at once. When a matrix of a new structure is encountered, the 340user calls spCreate(). This routine creates the basic frame for the 341linked-list and returns a pointer to this frame. This pointer is then 342passed as an argument to the other Sparse routines to indicate which matrix 343is to be operated on. The number of matrices that can be kept in memory at 344once is only limited by the amount of memory available to the user and the 345size of the matrices. When a matrix frame is no longer needed, the memory 346can be reclaimed by calling spDestroy(). 347 348 A more complete discussion of sparse systems of equations, methods for 349solving them, their error mechanisms, and the algorithms used in Sparse can 350be found in Kundert [kundert86]. A particular emphasis is placed on 351matrices resulting from circuit simulators. 352 353 3542.2: Error Control 355 356 There are two separate mechanisms that can cause errors during the 357factoring and solution of a system of equations. The first is ill- 358conditioning in the system. A system of equations is ill-conditioned if 359the solution is excessively sensitive to disturbances in the input data, 360which occurs when the system is nearly singular. If a system is ill- 361conditioned then uncertainty in the result is unavoidable, even if A is 362accurately factored into L and U. When ill-conditioning is a problem, the 363problem as stated is probably ill-posed and the system should be reformu- 364lated such that it is not so ill-conditioned. It is possible to measure 365the ill-conditioning of matrix using spCondition(). This function returns 366an estimate of the reciprocal of the condition number of the matrix (K(A)) 367[strang80]. The condition number can be used when computing a bound on the 368error in the solution using the following inequality [golub83]. 369 370 ||dx|| (||dA|| ||db||) 371 ------ < K(A) (------ + ------) + higher order terms 372 ||x|| (||A|| ||b|| ) 373 374where dA and db are the uncertainties in the matrix and right-hand side 375vector and are assumed small. 376 377 The second mechanism that causes uncertainty is the build up of round- 378off error. Roundoff error can become excessive if there is sufficient 379growth in the size of the elements during the factorization. Growth is 380controlled by careful pivoting. In Sparse, the pivoting is controlled by 381the relative threshold parameter. In conventional full matrix techniques 382the pivot is chosen to be the largest element in a column. When working 383with sparse matrices it is important to choose pivots to minimize the 384reduction in sparsity. The best pivot to retain sparsity is often not the 385best pivot to retain accuracy. Thus, some compromise must be made. In 386threshold pivoting, as used in this package, the best pivot to retain spar- 387sity is used unless it is smaller than the relative threshold times the 388largest element in the column. Thus, a relative threshold close to one 389emphasizes accuracy so it will produce a minimum amount of growth, unfor- 390tunately it also slows the factorization. A very small relative threshold 391emphasizes maintenance of sparsity and so speeds the factorization, but can 392result in a large amount of growth. In our experience, we have found that 393a relative threshold of 0.001 seems to result in a satisfactory compromise 394between speed and accuracy, though other authors suggest a more conserva- 395tive value of 0.1 [duff86]. 396 397 398 June 23, 1988 399 400 401 402 403 404 - 6 - 405 406 407 The growth that occurred during a factorization can be computed by 408taking the ratio of the largest matrix element in any stage of the factori- 409zation to the largest matrix element before factorization. The two numbers 410are estimated using spLargestElement(). If the growth is found to be 411excessive after spOrderAndFactor(), then the relative threshold should be 412increased and the matrix reconstructed and refactored. Once the matrix has 413been ordered and factored without suffering too much growth, the amount of 414growth that occurred should be recorded. If, on subsequent factorizations, 415as performed by spFactor(), the amount of growth becomes significantly 416larger, then the matrix should be reconstructed and reordered using the 417same relative threshold with spOrderAndFactor(). If the growth is still 418excessive, then the relative threshold should be raised again. 419 420 4212.3: Building the Matrix 422 423 It is not necessary to specify the size of the matrix before beginning 424to add elements to it. When the compiler option EXPANDABLE is turned on it 425is possible to initially specify the size of the matrix to any size equal 426to or smaller than the final size of the matrix. Specifically, the matrix 427size may be initially specified as zero. If this is done then, as the ele- 428ments are entered into the matrix, the matrix is enlarged as needed. This 429feature is particularly useful in circuit simulators because it allows the 430building of the matrix as the circuit description is parsed. Note that 431once the matrix has been reordered by the routines spMNA Preorder(), spFac- 432tor() or spOrderAndFactor() the size of the matrix becomes fixed and may no 433longer be enlarged unless the compiler option TRANSLATE is enabled. 434 435 The TRANSLATE option allows Sparse to translate a non-packed set of 436row and column numbers to an internal packed set. In other words, there 437may be rows and columns missing from the external description of the 438matrix. This feature provides two benefits. First, if two matrices are 439identical in structure, except for a few missing rows and columns in one, 440then the TRANSLATE option allows them to be treated identically. Simi- 441larly, rows and columns may be deleted from a matrix after it has been 442built and operated upon. Deletion of rows and columns is performed by the 443function spDeleteRowAndCol(). Second, it allows the use of the functions 444spGetElement(), spGetAdmittance(), spGetQuad(), and spGetOnes() after the 445matrix has been reordered. These functions access the matrix by using row 446and column indices, which have to be translated to internal indices once 447the matrix is reordered. Thus, when TRANSLATE is used in conjunction with 448the EXPANDABLE option, rows and columns may be added to a matrix after it 449has been reordered. 450 451 Another provided feature that is useful with circuit simulators is the 452ability to add elements to the matrix in row zero or column zero. These 453elements will have no affect on the matrix or the results. The benefit of 454this is that when working with a nodal formulation, grounded components do 455not have to be treated special when building the matrix. 456 457 458 459 460 461 462 463 464 June 23, 1988 465 466 467 468 469 470 - 7 - 471 472 4732.4: Initializing the Matrix 474 475 Once a matrix has been factored, it is necessary to clear the matrix 476before it can be reloaded with new values. The straight forward way of 477doing that is to call spClear(), which sets the value of every element in 478the matrix to zero. Sparse also provides a more flexible way to clear the 479matrix. Using spInitialize(), it is possible to clear and reload at least 480part of the matrix in one step. 481 482 Sparse allows the user to keep initialization information with each 483structurally nonzero matrix element. Each element has a pointer that is 484set and used by the user. The user can set this pointer using spInstallIn- 485itInfo() and may read it using spGetInitInfo(). The function spInitial- 486ize() is a user customizable way to initialize the matrix. Passed to this 487routine is a function pointer. spInitialize() sweeps through every element 488in the matrix and checks the pInitInfo pointer (the user supplied pointer). 489If the pInitInfo is NULL, which is true unless the user changes it (always 490true for fill-ins), then the element is zeroed. Otherwise, the function 491pointer is called and passed the pInitInfo pointer as well as the element 492pointer and the external row and column numbers, allowing the user to ini- 493tialize the matrix element and the right-hand side. 494 495 Why spInitialize() would be used over spClear() can be illustrated by 496way of an example. Consider a circuit simulator that handles linear and 497nonlinear resistors and capacitors performing a transient analysis. For 498the linear resistors, a constant value is loaded into the matrix at each 499time step and for each Newton iteration. For the linear capacitor, a value 500is loaded into the matrix that is constant over Newton iterations, but is a 501function of the time step and the integration method. The nonlinear com- 502ponents contribute values to the matrix that change on every time step and 503Newton iteration. 504 505 Sparse allows the user to attach a data structure to each element in 506the matrix. For this example, the user might attach a structure that held 507several pieces of information, such as the conductance of the linear resis- 508tor, the capacitance of the linear capacitor, the capacitance of the non- 509linear capacitor, and perhaps past values of capacitances. The user also 510provides a subroutine to spInitialize() that is called for each user- 511created element in the matrix. This routine would, using the information 512in the attached data structure, initialize the matrix element and perhaps 513the right-hand side vector. 514 515 In this example, the user supplied routine might load the linear con- 516ductance into the matrix and multiply it by some voltage to find a current 517that could be loaded into the right-hand side vector. For the capacitors, 518the routine would first apply an integration method and then load the 519matrix and the right-hand side. 520 521 This approach is useful for two reasons. First, much of the work of 522the device code in the simulator can be off-loaded onto the matrix package. 523Since there are usually many devices, this usually results overall in a 524simpler system. Second, the integration method can be hidden from the 525simulator device code. Thus the integration method can be changed simply 526by changing the routine handed to spInitialize(), resulting in a much 527 528 529 530 June 23, 1988 531 532 533 534 535 536 - 8 - 537 538 539cleaner and more easily maintained simulator. 540 541 5422.5: Indices 543 544 By far the most common errors made when using Sparse are related to 545array indices. Sparse itself contributes to the problem by having several 546different indexing schemes. There are three different options that affect 547index bounds or the way indices are interpreted. The first is 548ARRAY OFFSET, which only affects array indices. ARRAY OFFSET is a compiler 549flag that selects whether arrays start at index zero or index one. Note 550that if ARRAY OFFSET is zero then RHS[0] corresponds to row one in the 551matrix and Solution[0] corresponds to column one. Further note that when 552ARRAY OFFSET is set to one, then the allocated length of the arrays handed 553to the Sparse routines should be at least the external size of the matrix 554plus one. The main utility of ARRAY OFFSET is that it allows natural array 555indexing when Sparse is coupled to programs in other languages. For exam- 556ple; in FORTRAN arrays always start at one whereas in C array always start 557at zero. Thus the first entry in a FORTRAN array corresponds to the 558zero'th entry in a C array. Setting ARRAY OFFSET to zero allows the arrays 559in FORTRAN to start at one rather than two. For the rest of this discus- 560sion, assume that ARRAY OFFSET is set so that arrays start at one in the 561program that calls Sparse. 562 563 The second option that affects indices is EXPANDABLE. When EXPANDABLE 564is set false the upper bound on array and matrix indices is Size, where 565Size is a parameter handed to spCreate(). When EXPANDABLE set true, then 566there is essentially no upper bound on array indices. Indeed, the size of 567the matrix is determined by the largest row or column number handed to 568Sparse. The upper bound on the array indices then equals the final size 569determined by Sparse. This size can be determined by calling spGetSize(). 570 571 572 The final option that affects indices is TRANSLATE. This option was 573provided to allow row and columns to be deleted, but it also allows row and 574column numbers to be missing from the input description for a matrix. This 575means that the size of the matrix is not determined by the largest row or 576column number entered into the matrix. Rather, the size is determined by 577the total number of rows or column entered. For example, if the elements 578[2,3], [5,3], and [7,2] are entered into the matrix, the internal size of 579the matrix becomes four while the external size is seven. The internal 580size equals the number of rows and columns in the matrix while the external 581size equals the largest row or column number entered into the matrix. Note 582that if a row is entered into the matrix, then its corresponding column is 583also entered, and vice versa. The indices used in the RHS and Solution 584vectors correspond to the row and column indices in the matrix. Thus, for 585this example, valid data is expected in RHS at locations 2, 3, 5 and 7. 586Data at other locations is ignored. Similarly, valid data is returned in 587Solution at locations 2, 3, 5, and 7. The other locations are left 588unmolested. This shows that the length of the arrays correspond to the 589external size of the matrix. Again, this value can be determined by spGet- 590Size(). 591 592 593 594 595 596 June 23, 1988 597 598 599 600 601 602 - 9 - 603 604 6052.6: Configuring Sparse 606 607 It is possible at compile-time to customize Sparse for your particular 608application. This is done by changing the compiler options, which are kept 609in the personality file, spConfig.h. There are three classes of choices 610available. First are the Sparse options, which specify the dominant per- 611sonality characteristics, such as if real and/or complex systems of equa- 612tions are to be handled. The second class is the Sparse constants, such as 613the default pivot threshold and the amount of memory initially allocated 614per matrix. The last class is the machine constants. These numbers must 615be updated when Sparse is ported to another machine. 616 617 As an aid in the setup and testing of Sparse a test routine and 618several test matrices and their solutions have been provided. The test 619routine is capable of reading files generated by spFileMatrix() and 620spFileVector(). 621 622 By default Sparse stores all real numbers and performs all computa- 623tions using double precision arithmetic. This can be changed by changing 624the definition of spREAL from double to float. spREAL is defined in 625spExports.h. 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 June 23, 1988 663 664 665 666 667 668 - 10 - 669 670 6713: INTRODUCTION TO THE SPARSE ROUTINES 672 673In this section the routines are grouped by function and briefly described. 674 6753.1: Creating the Matrix 676 677spCreate() 678 Allocates and initializes the data structure for a matrix. Necessari- 679 ly the first routine run for any particular matrix. 680 681spDestroy() 682 Destroys the data structure for a matrix and frees the memory. 683 684spSetReal() 685spSetComplex() 686 These routines toggle a flag internal to Sparse that indicates that 687 the matrix is either real or complex. This is useful if both real and 688 complex matrices of identical structure are expected. 689 690 6913.2: Building the Matrix 692 693spGetElement() 694 Assures that the specified element exists in the matrix data structure 695 and returns a pointer to it. 696 697spGetAdmittance() 698spGetQuad() 699spGetOnes() 700 These routines add a group of four related elements to the matrix. 701 spGetAdmittance() adds the four elements associated with a two termi- 702 nal admittance. spGetQuad() is a more general routine that is useful 703 for entering controlled sources to the matrix. spGetOnes() adds the 704 four structural ones to the matrix that are often encountered with 705 elements that do not have admittance representations. 706 707spDeleteRowAndCol() 708 This function is used to delete a row and column from the matrix. 709 710 7113.3: Clearing the Matrix 712 713spClear() 714 Sets every element in the matrix to zero. 715 716spInitialize() 717 Runs a user provided initialization routine on each element in the ma- 718 trix. This routine would be used in lieu of spClear(). 719 720spGetInitInfo() 721spInstallInitInfo() 722 These routines allow the user to install and read a user-provided 723 pointer to initialization data for a particular matrix element. 724 725 726 727 728 June 23, 1988 729 730 731 732 733 734 - 11 - 735 736 737 738spStripFills() 739 This routine returns a matrix to a semi-virgin state by removing all 740 fill-ins. This can be useful if a matrix is to be reordered and it 741 has changed significantly since it was previously ordered. This may 742 be the case if a few rows and columns have been added or deleted or if 743 the previous ordering was done on a matrix that was numerically quite 744 different than the matrix currently being factored. Stripping and 745 reordering a matrix may speed subsequent factorization if the current 746 ordering is inferior, whereas simply reordering will generally only 747 enhance accuracy and not speed. 748 749 7503.4: Placing Data in the Matrix 751 752spADD REAL ELEMENT() 753spADD IMAG ELEMENT() 754spADD COMPLEX ELEMENT() 755 Adds a value to a particular matrix element. 756 757spADD REAL QUAD() 758spADD IMAG QUAD() 759spADD COMPLEX QUAD() 760 Adds a value to a group of four matrix elements. 761 762 7633.5: Influencing the Factorization 764 765spMNA Preorder() 766 This routine preorders modified node admittance matrices so that 767 Sparse can take full advantage of their structure. In particular, 768 this routine tries to remove zeros from the diagonal so that diagonal 769 pivoting can be used more successfully. 770 771spPartition() 772 Sparse partitions the matrix in an attempt to make spFactor() run as 773 fast as possible. The partitioning is a relatively expensive opera- 774 tion that is not needed in all cases. spPartition() allows the user 775 specify a simpler and faster partitioning. 776 777spScale() 778 It is sometimes desirable to scale the rows and columns of a matrix in 779 to achieve a better pivoting order. This is particularly true in 780 modified node admittance matrices, where the size of the elements in a 781 matrix can easily vary through ten to twelve orders of magnitude. 782 This routine performs scaling on a matrix. 783 784 7853.6: Factoring the Matrix 786 787spOrderAndFactor() 788 This routine chooses a pivot order for the matrix and factors it into 789 LU form. It handles both the initial factorization and subsequent 790 factorizations when a reordering is desired. 791 792 793 794 June 23, 1988 795 796 797 798 799 800 - 12 - 801 802 803 804spFactor() 805 Factors a matrix that has already been ordered by spOrderAndFactor(). 806 If spFactor() is passed a matrix that needs ordering, it will automat- 807 ically pass the matrix to spOrderAndFactor(). 808 809 8103.7: Solving the Matrix Equation 811 812spSolve() 813 Solves the matrix equation 814 815 Ax = b 816 given the matrix A factored into LU form and b. 817 818spSolveTransposed() 819 When working with adjoint systems, such as in sensitivity analysis, it 820 is desirable to quickly solve 821 822 T 823 A x = b 824 Once A has been factored into LU form, this routine can be used to 825 solve the transposed system without having to suffer the cost of fac- 826 toring the matrix again. 827 828 8293.8: Numerical Error Estimation 830 831spCondition() 832 Estimates the L-infinity condition number of the matrix. This number 833 is a measure of the ill-conditioning in the matrix equation. It is 834 also useful for making estimates of the error in the solution. 835 836spNorm() 837 Returns the L-infinity norm (the maximum absolute row sum) of an un- 838 factored matrix. 839 840spPseudoCondition() 841 Returns the ratio of the largest pivot to the smallest pivot of a fac- 842 tored matrix. This is a rough indicator of ill-conditioning in the 843 matrix. 844 845spLargestElement() 846 If passed an unfactored matrix, this routine returns the absolute 847 value of the largest element in the matrix. If passed a factored ma- 848 trix, it returns an estimate of the largest element that occurred in 849 any of the reduced submatrices during the factorization. The ratio of 850 these two numbers (factored/unfactored) is the growth, which is used 851 to determine if the pivoting order is numerically satisfactory. 852 853spRoundoff() 854 Returns a bound on the magnitude of the largest element in E = A-LU, 855 where E represents error in the matrix resulting from roundoff error 856 during the factorization. 857 858 859 860 861 June 23, 1988 862 863 864 865 866 867 - 13 - 868 869 8703.9: Matrix Operations 871 872spDeterminant() 873 This routine simply calculates and returns the determinant of the fac- 874 tored matrix. 875 876spMultiply() 877 This routine multiplys the matrix by a vector on the right. This is 878 useful for forming the product Ax = b in order to determine if a cal- 879 culated solution is correct. 880 881spMultTransposed() 882 Multiplys the transposed matrix by a vector on the right. This is 883 useful for forming the product A sup {roman T} x = b in order to 884 determine if a calculated solution is correct. 885 886 8873.10: Matrix Statistics and Documentation 888 889spError() 890 Determines the error status of a particular matrix. While most of the 891 Sparse routines do return an indication that an error has occurred, 892 some do not and so spError() provides the only way of uncovering these 893 errors. 894 895spWhereSingular() 896 Returns the row and column number where the matrix was detected as 897 singular or where a zero pivot was found. 898 899spGetSize() 900 A function that returns the size of the matrix. Either the internal 901 or external size of the matrix is returned. The internal size of the 902 matrix is the actual size of the matrix whereas the external size is 903 the value of the largest row or column number. These two numbers may 904 differ if the TRANSLATE option is used. 905 906spElementCount() 907spFillinCount() 908 Functions that return the total number of elements in the matrix, and 909 the number of fill-ins in the matrix. These functions are useful for 910 gathering statistics on matrices. 911 912spPrint() 913 This routine outputs the matrix as well as some statistics to standard 914 output in a format that is readable by people. The matrix can be 915 printed in either a compressed or standard format. In the standard 916 format, a numeric value is given for each structurally nonzero ele- 917 ment, whereas in the compressed format, only the existence or nonex- 918 istence of an element is indicated. This routine is not suitable for 919 use on large matrices. 920 921 922 923 924 925 926 927 June 23, 1988 928 929 930 931 932 933 - 14 - 934 935 936 937spFileMatrix() 938spFileVector() 939 These two routines send a copy of the matrix and its right-hand side 940 vector to a file. This file can then be read by the test program that 941 is included with Sparse. Only those elements of the matrix that are 942 structurally nonzero are output, so very large matrices can be sent to 943 a file. 944 945spFileStats() 946 This routine calculates and sends some useful statistics concerning a 947 matrix to a file. 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 June 23, 1988 994 995 996 997 998 999 - 15 - 1000 1001 10024: SPARSE ROUTINES 1003 1004This section contains a complete list of the Sparse routines that are 1005available to the user. Each routine is described as to its function and 1006how to use it. The routines are listed in alphabetic order. 1007 1008 1009 1010 1011 10124.1: spClear() 1013 1014Sets every element in the matrix to zero. The Sparse error state is 1015cleared to spOKAY in this routine. 1016 1017void spClear( Matrix ) 1018 1019o Argument: 1020 1021 Matrix input (char *) 1022 Pointer to matrix that is to be cleared. 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 June 23, 1988 1061 1062 1063 1064 1065 1066 - 16 - 1067 1068 1069 1070 1071 10724.2: spCondition() 1073 1074spCondition() computes an estimate of the condition number using a varia- 1075tion on the LINPACK condition number estimation algorithm. This quantity 1076is an measure of ill-conditioning in the matrix. To avoid problems with 1077overflow, the reciprocal of the condition number is returned. If this 1078number is small, and if the matrix is scaled such that uncertainties in the 1079RHS and the matrix entries are equilibrated, then the matrix is ill- 1080conditioned. If the this number is near one, the matrix is well condi- 1081tioned. This routine must only be used after a matrix has been factored by 1082spOrderAndFactor() or spFactor() and before it is cleared by spClear() or 1083spInitialize(). 1084 1085Unlike the LINPACK condition number estimator, this routines returns the L 1086infinity condition number. This is an artifact of Sparse placing ones on 1087the diagonal of the upper triangular matrix rather than the lower. This 1088difference should be of no importance. 1089 1090spREAL spCondition( Matrix, NormOfMatrix, Error ) 1091 1092o Returns: 1093 An estimate of the L infinity condition number of the matrix. 1094 1095o Arguments: 1096 1097 Matrix input (char *) 1098 The matrix for which the condition number is desired. 1099 1100 NormOfMatrix input (spREAL) 1101 The L-infinity norm of the unfactored matrix as computed by 1102 spNorm(). 1103 1104 Error output (int *) 1105 Used to return the error code. 1106 1107o Possible errors: 1108 spSINGULAR 1109 spNO MEMORY 1110 Error is not cleared in this routine. 1111 1112o Compiler options that must be set for this routine to exist: 1113 CONDITION 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 June 23, 1988 1128 1129 1130 1131 1132 1133 - 17 - 1134 1135 1136 1137 1138 11394.3: spCreate() 1140 1141Allocates and initializes the data structures associated with a matrix. 1142This routine is necessarily the first routine run for any particular ma- 1143trix. 1144 1145char *spCreate( Size, Complex, Error ) 1146 1147o Returned: 1148 A pointer to the matrix is returned cast into the form of a pointer to 1149 a character. This pointer is then passed and used by the other matrix 1150 routines to refer to a particular matrix. If an error occurs, the 1151 NULL pointer is returned. 1152 1153o Arguments: 1154 1155 Size input (int) 1156 Size of matrix. When the compiler option EXPANDABLE is turned 1157 on, Size is used as a lower bound on the size of the matrix. 1158 Size must not be negative. 1159 1160 Complex input (int) 1161 Type of matrix. If Complex is 0 then the matrix is real, other- 1162 wise the matrix will be complex. Note that if the routines are 1163 not set up to handle the type of matrix requested, then a spPANIC 1164 error will occur. 1165 1166 Error output (int *) 1167 Returns error flag, needed because function spError() will not 1168 work correctly if spCreate() returns NULL. 1169 1170o Possible errors: 1171 spNO MEMORY 1172 spPANIC 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 June 23, 1988 1195 1196 1197 1198 1199 1200 - 18 - 1201 1202 1203 1204 1205 12064.4: spDeleteRowAndCol() 1207 1208This function is used to delete a row and column from the matrix. The ele- 1209ments removed from the matrix are never used again and are not freed until 1210the matrix is destroyed and so the pointers to these elements remain valid. 1211 1212void spDeleteRowAndCol( Matrix, Row, Col ) 1213 1214o Arguments: 1215 1216 Matrix input (char *) 1217 The matrix from which the row and column are to be deleted. 1218 1219 Row input (int) 1220 The row to be deleted. 1221 1222 Col input (int) 1223 The column to be deleted. 1224 1225o Compiler options that must be set for this routine to exist: 1226 DELETE 1227 TRANSLATE 1228 1229 1230 1231 1232 12334.5: spDestroy() 1234 1235Destroys a matrix frame and reclaims the memory. 1236 1237void spDestroy( Matrix ) 1238 1239o Argument: 1240 1241 Matrix input (char *) 1242 Pointer to the matrix frame which is to be removed from memory. 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 June 23, 1988 1263 1264 1265 1266 1267 1268 - 19 - 1269 1270 1271 1272 1273 12744.6: spDeterminant() 1275 1276This routine in capable of calculating the determinant of the matrix once 1277the LU factorization has been performed. Hence, only use this routine 1278after spFactor() or spOrderAndFactor() and before spClear() or spInitial- 1279ize(). Note that the determinants of matrices can be very large or very 1280small. On large matrices, the determinant can be far larger or smaller 1281than can be represented by a floating point number. For this reason the 1282mantissa and exponent of the determinant are returned separately. 1283 1284void spDeterminant( Matrix, Exponent, Determinant ) 1285void spDeterminant( Matrix, Exponent, Determinant, iDeterminant ) 1286 1287o Arguments: 1288 1289 Matrix input (char *) 1290 The matrix for which the determinant is desired. 1291 1292 Exponent output (int *) 1293 The logarithm base 10 of the scale factor for the determinant. 1294 To find the actual determinant, Exponent should be added to the 1295 exponent of Determinant and iDeterminant. 1296 1297 Determinant output (spREAL *) 1298 The real portion of the determinant. If the matrix is real, then 1299 the magnitude of this number is scaled to be greater than or 1300 equal to 1.0 and less than 10.0. Otherwise the magnitude of the 1301 complex determinant will be scaled such. 1302 1303 iDeterminant output (spREAL *) 1304 The imaginary portion of the determinant. When the matrix is 1305 real this pointer need not be supplied; nothing will be returned. 1306 1307o Compiler options that must be set for this routine to exist: 1308 DETERMINANT 1309 1310o Bugs: 1311 The sign of determinant may be in error if rows and columns have been 1312 added or deleted from matrix. 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 June 23, 1988 1330 1331 1332 1333 1334 1335 - 20 - 1336 1337 1338 1339 1340 13414.7: spElementCount() 1342 1343Returns the total number of structurally nonzero elements in the matrix. 1344 1345int spElementCount( Matrix ) 1346 1347o Returns: 1348 The total number of structurally nonzero elements. 1349 1350o Argument: 1351 1352 Matrix input (char *) 1353 Pointer to the matrix. 1354 1355 1356 1357 1358 13594.8: spError() 1360 1361This function returns the error status of a matrix. 1362 1363int MatrixError( Matrix ) 1364 1365o Returned: 1366 The error status of the given matrix. 1367 1368o Argument: 1369 1370 Matrix input (char *) 1371 The matrix for which the error status is desired. 1372 1373o Possible errors: 1374 spOKAY 1375 spILL CONDITIONED 1376 spZERO PIVOT 1377 spSINGULAR 1378 spNO MEMORY 1379 spPANIC 1380 Error is not cleared in this routine. 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 June 23, 1988 1398 1399 1400 1401 1402 1403 - 21 - 1404 1405 1406 1407 1408 14094.9: spFactor() 1410 1411This routine factors the matrix into LU form and is the companion routine 1412to spOrderAndFactor(). Unlike spOrderAndFactor(), spFactor() cannot change 1413the ordering. Its utility is that it is considerably faster. The standard 1414way to use these two routines is to first use spOrderAndFactor() for the 1415initial factorization. For subsequent factorizations, spFactor() is used. 1416If spFactor() is called for the initial factorization of the matrix, then 1417it will automatically call spOrderAndFactor() with the default thresholds. 1418If spFactor() finds a zero on the diagonal, it will terminate early and 1419complain. This does not necessarily mean that matrix is singular. Before 1420a matrix is condemned as being singular, it should be run through spOr- 1421derAndFactor(), which can reorder the matrix and remove the offensive zero 1422from the diagonal. 1423 1424int spFactor( Matrix ) 1425 1426o Returned: 1427 The error code is returned. Possible errors are listed below. 1428 1429o Argument: 1430 1431 Matrix input (char *) 1432 Pointer to matrix to be factored. 1433 1434o Possible errors: 1435 spZERO PIVOT 1436 spNO MEMORY 1437 spSINGULAR 1438 spILL CONDITIONED 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 June 23, 1988 1465 1466 1467 1468 1469 1470 - 22 - 1471 1472 1473 1474 14754.10: spFileMatrix() 1476 1477Writes matrix to file in format suitable to be read back in by the matrix 1478test program. Normally, spFileMatrix() should be executed before the ma- 1479trix is factored, otherwise matrix is output in factored form. If the ma- 1480trix is sent to a file without the header or data, it will be in a form 1481that is easily plotted by typical plotting programs. 1482 1483int spFileMatrix( Matrix, File, Label, Reordered, Data, Header ) 1484 1485o Returns: 1486 One is returned if routine was successful, otherwise zero is returned. 1487 The calling function can query errno (the system global error vari- 1488 able) as to the reason why this routine failed. 1489 1490o Arguments: 1491 1492 Matrix input (char *) 1493 Pointer to matrix that is to be sent to file. 1494 1495 File input (char *) 1496 Name of output file. 1497 1498 Label input (char *) 1499 String that is transferred to file and used as a label. String 1500 should fit on one line and have no embedded line feeds. 1501 1502 Reordered input (int) 1503 Specifies whether the matrix should be output using the original 1504 order or in reordered form. Zero specifies original order. 1505 1506 Data input (int) 1507 Indicates that the element values should be output along with the 1508 indices for each element. Element values are not output if Data 1509 is zero. This parameter must be nonzero if matrix is to be read 1510 by the Sparse test program. 1511 1512 Header input (int) 1513 If nonzero a header is output that includes that size of the ma- 1514 trix and the label. This parameter must be nonzero if matrix is 1515 to be read by the Sparse test program. 1516 1517o Compiler options that must be set for this routine to exist: 1518 DOCUMENTATION 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 June 23, 1988 1531 1532 1533 1534 1535 1536 - 23 - 1537 1538 1539 1540 1541 15424.11: spFileStats() 1543 1544Appends useful information concerning the matrix to the end of a file. If 1545file does not exist, it is created. This file should not be the same as 1546one used to hold the matrix or vector if the matrix is to be read by the 1547Sparse test program. Should be executed after the matrix is factored. 1548 1549int spFileStats( Matrix, File, Label ) 1550 1551o Returns: 1552 One is returned if routine was successful, otherwise zero is returned. 1553 The calling function can query errno (the system global error vari- 1554 able) as to the reason why this routine failed. 1555 1556o Arguments: 1557 1558 Matrix input (char *) 1559 Pointer to matrix for which statistics are desired. 1560 1561 File input (char *) 1562 Name of output file. 1563 1564 Label input (char *) 1565 String that is transferred to file and is used as a label. String 1566 should fit on one line and have no embedded line feeds. 1567 1568o Compiler options that must be set for this routine to exist: 1569 DOCUMENTATION 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 June 23, 1988 1598 1599 1600 1601 1602 1603 - 24 - 1604 1605 1606 1607 1608 16094.12: spFileVector() 1610 1611Appends the RHS vector to the end of a file in a format suitable to be read 1612back in by the matrix test program. If file does not exist, it is created. 1613To be compatible with the test program, if spFileVector() is run, it must 1614be run after spFileMatrix() and use the same file. 1615 1616int spFileVector( Matrix, File, RHS ) 1617int spFileVector( Matrix, File, RHS, iRHS ) 1618 1619o Returns: 1620 One is returned if routine was successful, otherwise zero is returned. 1621 The calling function can query errno (the system global error vari- 1622 able) as to the reason why this routine failed. 1623 1624o Arguments: 1625 1626 Matrix input (char *) 1627 Pointer to matrix that corresponds to the vector to be output. 1628 1629 File input (char *) 1630 Name of file where output is to be written. 1631 1632 RHS input (spREAL[]) 1633 The right-hand side vector. RHS contains only the real portion 1634 of the right-hand side vector if the matrix is complex and 1635 spSEPARATED COMPLEX VECTORS is set true. 1636 1637 iRHS input (spREAL[]) 1638 Right-hand side vector, imaginary portion. Not necessary if ma- 1639 trix is real or if spSEPARATED COMPLEX VECTORS is set false. 1640 1641o Compiler options that must be set for this routine to exist: 1642 DOCUMENTATION 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 June 23, 1988 1665 1666 1667 1668 1669 1670 - 25 - 1671 1672 1673 1674 1675 16764.13: spFillinCount() 1677 1678Returns the total number of fill-ins in the matrix. A fill-in is an ele- 1679ment that is originally structurally zero, but becomes nonzero during the 1680factorization. 1681 1682int spFillinCount( Matrix ) 1683 1684o Returns: 1685 The total number of fill-ins. 1686 1687o Argument: 1688 1689 Matrix input (char *) 1690 Pointer to the matrix. 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 June 23, 1988 1732 1733 1734 1735 1736 1737 - 26 - 1738 1739 1740 1741 1742 17434.14: spGetAdmittance() 1744 1745Performs same function as spGetElement() except rather than one element, 1746all four matrix elements for a floating admittance are reserved. This rou- 1747tine also works if the admittance is grounded (zero is the ground node). 1748This function returns a group of pointers to the four elements through Tem- 1749plate, which is an output. They are used by the spADD QUAD() macros to 1750directly access matrix elements during subsequent loads of the matrix. 1751spGetAdmittance() arranges the pointers in Template so that the 1752spADD QUAD() routines add the admittance to the elements at [Node1,Node1] 1753and [Node2,Node2] and subtract the admittance from the elements at 1754[Node1,Node2] and [Node2,Node1]. This routine is only to be used before 1755spMNA Preorder(), spFactor() or spOrderAndFactor() unless the compiler flag 1756TRANSLATE is enabled. 1757 1758int spGetAdmittance( Matrix, Node1, Node2, Template ) 1759 1760o Returned: 1761 The error code is returned. Possible errors are listed below. 1762 spGetAdmittance() does not clear the error state, so it is possible to 1763 ignore the return code of each spGetAdmittance() call, and check for 1764 errors after constructing the whole matrix by calling spError(). 1765 1766o Arguments: 1767 1768 Matrix input (char *) 1769 Pointer to the matrix that admittance is to be installed. 1770 1771 Node1 input (int) 1772 One node number for the admittance. Node1 must be in the range 1773 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 1774 flags are set true. In either case Node1 must not be negative. 1775 1776 Node2 input (int) 1777 Other node number for the admittance. Node2 must be in the range 1778 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 1779 flags are set true. In either case Node2 must not be negative. 1780 1781 Template output (struct spTemplate *) 1782 Collection of pointers to four elements that are later used to 1783 directly address elements. User must supply the template, this 1784 routine will fill it. 1785 1786o Possible errors: 1787 spNO MEMORY 1788 Error is not cleared in this routine. 1789 1790o Compiler options that must be set for this routine to exist: 1791 QUAD ELEMENT 1792 1793 1794 1795 1796 1797 1798 June 23, 1988 1799 1800 1801 1802 1803 1804 - 27 - 1805 1806 1807 1808 1809 18104.15: spGetElement() 1811 1812Reserves an element at [Row,Col] and returns a pointer to it. If element 1813is not found then it is created and spliced into matrix. A pointer to the 1814real portion of the element is returned. This pointer is later used by the 1815spADD ELEMENT() macros to directly access the element. This routine is 1816only to be used before spMNA Preorder(), spFactor() or spOrderAndFactor() 1817unless the compiler option TRANSLATE is set true. 1818 1819spREAL *spGetElement( Matrix, Row, Col ) 1820 1821o Returned: 1822 Returns a pointer to the element. This pointer is then used to 1823 directly access the element during successive builds. Returns NULL if 1824 insufficient memory is available. spGetElement() does not clear the 1825 error state, so it is possible to ignore the return code of each 1826 spGetElement() call, and check for errors after constructing the whole 1827 matrix by calling spError(). 1828 1829o Arguments: 1830 1831 Matrix input (char *) 1832 Pointer to the matrix that the element is to be added to. 1833 1834 Row input (int) 1835 Row index for element. Row must be in the range [0..Size] unless 1836 either the TRANSLATE or EXPANDABLE compiler flags are set true. 1837 In either case Row must not be negative though it may be zero. 1838 If zero then the element is not entered into the matrix, but is 1839 otherwise treated normally. 1840 1841 Col input (int) 1842 Column index for element. Col must be in the range [0..Size] un- 1843 less either the TRANSLATE or EXPANDABLE compiler flags are set 1844 true. In either case Col must not be negative though it may be 1845 zero. If zero then the element is not entered into the matrix, 1846 but is otherwise treated normally. 1847 1848o Possible errors: 1849 spNO MEMORY 1850 Error is not cleared in this routine. 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 June 23, 1988 1866 1867 1868 1869 1870 1871 - 28 - 1872 1873 1874 1875 1876 18774.16: spGetInitInfo() 1878 1879With the INITIALIZE compiler option enabled Sparse allows the user to keep 1880initialization information with each structurally nonzero matrix element. 1881Each element has a pointer (referred to as pInitInfo) that is set and used 1882by the user. This routine returns pInitInfo from a particular matrix ele- 1883ment. 1884 1885char *spGetInitInfo( pElement ) 1886 1887o Returned: 1888 The user installed pointer pInitInfo. 1889 1890o Argument: 1891 1892 pElement input (spREAL *) 1893 Pointer to the element to which pInitInfo is attached. 1894 1895o Compiler options that must be set for this routine to exist: 1896 INITIALIZE 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 June 23, 1988 1933 1934 1935 1936 1937 1938 - 29 - 1939 1940 1941 1942 19434.17: spGetOnes() 1944 1945Performs a similar function to spGetAdmittance() except that the four 1946reserved matrix elements are assumed to be structural ones generated by 1947components without admittance representations during a modified nodal 1948analysis. Positive ones are placed at [Pos,Eqn] and [Eqn,Pos] and negative 1949ones are placed at [Neg,Eqn] and [Eqn,Neg]. This function returns a group 1950of pointers to the four elements through Template, which is an output. 1951They are used by the spADD QUAD() macros to add the ones directly to the 1952matrix elements during subsequent loads of the matrix. This routine is 1953only to be used before spMNA Preorder(), spFactor() or spOrderAndFactor() 1954unless the compiler flag TRANSLATE is set true. 1955 1956int spGetOnes( Matrix, Pos, Neg, Eqn, Template ) 1957 1958o Returned: 1959 The error code is returned. Possible errors are listed below. 1960 spGetOnes() does not clear the error state, so it is possible to ig- 1961 nore the return code of each spGetOnes() call, and check for errors 1962 after constructing the whole matrix by calling spError(). 1963 1964o Arguments: 1965 1966 Matrix input (char *) 1967 Pointer to the matrix that ones are to be entered in. 1968 1969 Pos input (int) 1970 Number of positive node. Must be in the range of [0..Size] un- 1971 less either the options EXPANDABLE or TRANSLATE are used. Zero 1972 is the ground row. In no case may Pos be less than zero. 1973 1974 Neg input (int) 1975 Number of negative node. Must be in the range of [0..Size] un- 1976 less either the options EXPANDABLE or TRANSLATE are used. Zero is 1977 the ground row. In no case may Neg be less than zero. 1978 1979 Eqn input (int) 1980 Row that contains the branch equation. Must be in the range of 1981 [1..Size] unless either the options EXPANDABLE or TRANSLATE are 1982 used. In no case may Eqn be less than one. 1983 1984 Template output (struct spTemplate *) 1985 Collection of pointers to four elements that are later used to 1986 directly address elements. User must supply the template, this 1987 routine will fill it. 1988 1989o Possible errors: 1990 spNO MEMORY 1991 Error is not cleared in this routine. 1992 1993o Compiler options that must be set for this routine to exist: 1994 QUAD ELEMENT 1995 1996 1997 1998 June 23, 1988 1999 2000 2001 2002 2003 2004 - 30 - 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 June 23, 1988 2065 2066 2067 2068 2069 2070 - 31 - 2071 2072 2073 2074 2075 20764.18: spGetQuad() 2077 2078Similar to spGetAdmittance(), except that spGetAdmittance() only handles 20792-terminal components, whereas spGetQuad() handles simple 4-terminals as 2080well. These 4-terminals are simply generalized 2-terminals with the option 2081of having the sense terminals different from the source and sink terminals. 2082spGetQuad() installs four elements into the matrix and returns their 2083pointers in the Template structure, which is an output. The pointers are 2084arranged in Template such that when passed to one of the spADD QUAD() mac- 2085ros along with an admittance, the admittance will be added to the elements 2086at [Row1,Col1] and [Row2,Col2] and subtracted from the elements at 2087[Row1,Col2] and [Row2,Col1]. The routine works fine if any of the rows and 2088columns are zero. This routine is only to be used before spMNA Preorder(), 2089spFactor() or spOrderAndFactor() unless TRANSLATE is set true. 2090 2091int spGetQuad( Matrix, Row1, Row2, Col1, Col2, Template ) 2092 2093o Returned: 2094 The error code is returned. Possible errors are listed below. spGet- 2095 Quad() does not clear the error state, so it is possible to ignore the 2096 return code of each spGetQuad() call, and check for errors after con- 2097 structing the whole matrix by calling spError(). 2098 2099o Arguments: 2100 2101 Matrix input (char *) 2102 Pointer to the matrix that quad is to be entered in. 2103 2104 Row1 input (int) 2105 First row index for the elements. Row1 must be in the range 2106 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 2107 flags are set true. In either case Row1 must not be negative. 2108 2109 Row2 input (int) 2110 Second row index for the elements. Row2 must be in the range 2111 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 2112 flags are set true. In either case Row2 must not be negative. 2113 2114 Col1 input (int) 2115 First column index for the elements. Col1 must be in the range 2116 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 2117 flags are set true. In either case Col1 must not be negative. 2118 2119 Col2 input (int) 2120 Second column index for the elements. Col2 must be in the range 2121 [0..Size] unless either the TRANSLATE or EXPANDABLE compiler 2122 flags are set true. In either case Col2 must not be negative. 2123 2124 Template output (struct spTemplate *) 2125 Collection of pointers to four elements that are later used to 2126 directly address elements. User must supply the template, this 2127 routine will fill it. 2128 2129 2130 2131 June 23, 1988 2132 2133 2134 2135 2136 2137 - 32 - 2138 2139 2140o Possible errors: 2141 spNO MEMORY 2142 Error is not cleared in this routine. 2143 2144o Compiler options that must be set for this routine to exist: 2145 QUAD ELEMENT 2146 2147 2148 2149 2150 21514.19: spGetSize() 2152 2153Returns the size of the matrix, either the internal or external size of the 2154matrix is returned. The internal size is the actual number of rows and 2155columns in the matrix. The external size is equal to the largest row or 2156column number. These numbers will be the same unless the TRANSLATE option 2157is enabled. 2158 2159int spGetSize( Matrix, External ) 2160 2161o Returned: 2162 The size of the matrix. 2163 2164o Arguments: 2165 2166 Matrix input (char *) 2167 Pointer to the matrix for which the size is desired. 2168 2169 External input (int) 2170 If External is nonzero, the external size of the matrix is re- 2171 turned, otherwise the internal size of the matrix is returned. 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 June 23, 1988 2199 2200 2201 2202 2203 2204 - 33 - 2205 2206 2207 2208 2209 22104.20: spInitialize() 2211 2212spInitialize() is a user customizable way to initialize the matrix. Passed 2213to this routine is a function pointer. spInitialize() sweeps through every 2214element in the matrix and checks the pInitInfo pointer (the user supplied 2215pointer). If the pInitInfo is NULL, which is true unless the user changes 2216it (always true for fill-ins), then the element is zeroed. Otherwise, the 2217function pointer is called and passed the pInitInfo pointer as well as the 2218element pointer and the external row and column numbers allowing the user 2219to set the value of each element and perhaps the right-hand side vector. 2220 2221The user function (pInit()) is expected to return a nonzero integer if 2222there is a fatal error and zero otherwise. Upon encountering a nonzero re- 2223turn code, spInitialize() terminates and returns the error code. 2224 2225The Sparse error state is cleared to spOKAY in this routine. 2226 2227int spInitialize( Matrix, pInit ) 2228 2229o Returns: 2230 The error code returned by pInit. 2231 2232o Arguments: 2233 2234 Matrix input (char *) 2235 Pointer to the matrix that is to be initialized. 2236 2237 pInit input ((*int)()) 2238 Pointer to a function that, given a pointer to an element, a 2239 pointer to the users data structure containing initialization in- 2240 formation for that element, and the row and column number of the 2241 element, initializes it. 2242 2243 2244int pInit( pElement, pInitInfo, Row, Col ) 2245 2246o Returns: 2247 Nonzero if fatal error, zero otherwise. 2248 2249o Arguments: 2250 2251 pElement input (spREAL *) 2252 The pointer to the real portion of the element. The real portion 2253 can be accessed using either *pElement or pElement[0]. The ima- 2254 ginary portion can be accessed using either *(pElement+1) or 2255 pElement[1]. 2256 2257 pInitInfo input (char *) 2258 The user-installed pointer to the initialization data structure. 2259 2260 Row input (int) 2261 The external row number of the element. 2262 2263 2264 2265 June 23, 1988 2266 2267 2268 2269 2270 2271 - 34 - 2272 2273 2274 Col input (int) 2275 The external column number of the element. 2276 2277o Compiler options that must be set for this routine to exist: 2278 INITIALIZE 2279 2280 2281 2282 2283 22844.21: spInstallInitInfo() 2285 2286With the INITIALIZE compiler option enabled Sparse allows the user to keep 2287initialization information with each structurally nonzero matrix element. 2288Each element has a pointer (referred to as pInitInfo) that is set and used 2289by the user. This routine installs the pointer pInitInfo into a particular 2290matrix element. 2291 2292void spInstallInitInfo( pElement, pInitInfo ) 2293 2294o Arguments: 2295 2296 pElement input (spREAL *) 2297 Pointer to the element to which pInitInfo is to be attached. 2298 2299 pInitInfo input (char *) 2300 The pointer pInitInfo. 2301 2302o Compiler options that must be set for this routine to exist: 2303 INITIALIZE 2304 2305 2306 2307 2308 2309 2310 2311 2312 2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 2325 2326 2327 2328 2329 2330 2331 2332 June 23, 1988 2333 2334 2335 2336 2337 2338 - 35 - 2339 2340 2341 2342 2343 23444.22: spLargestElement() 2345 2346If this routine is called before the matrix is factored, it returns the ab- 2347solute value of the largest element in the matrix. If called after the ma- 2348trix has been factored, it returns a lower bound on the absolute value of 2349the largest element that occurred in any of the reduced submatrices during 2350the factorization. The ratio of these two numbers (factored/unfactored) is 2351the growth, which can be used to determine if the pivoting order is ade- 2352quate. A large growth implies that considerable error has been made in the 2353factorization and that it is probably a good idea to reorder the matrix. 2354If a large growth in encountered after using spFactor(), reconstruct the 2355matrix and refactor using spOrderAndFactor(). If a large growth is encoun- 2356tered after using spOrderAndFactor(), refactor using spOrderAndFactor() 2357with the pivot threshold increased, say to 0.1. 2358 2359spREAL spLargestElement( Matrix ) 2360 2361o Returns: 2362 If matrix is unfactored, returns the magnitude of the largest element 2363 in the matrix. If the matrix is factored, a bound on the magnitude of 2364 the largest element in any of the reduced submatrices is returned. 2365 2366o Argument: 2367 2368 Matrix input (char *) 2369 Pointer to the matrix. 2370 2371o Compiler options that must be set for this routine to exist: 2372 STABILITY 2373 2374 2375 2376 2377 2378 2379 2380 2381 2382 2383 2384 2385 2386 2387 2388 2389 2390 2391 2392 2393 2394 2395 2396 2397 2398 2399 June 23, 1988 2400 2401 2402 2403 2404 2405 - 36 - 2406 2407 2408 2409 2410 24114.23: spMNA Preorder() 2412 2413This routine massages modified node admittance matrices to improve the per- 2414formance of spOrderAndFactor(). It tries to remove structural zeros from 2415the diagonal by exploiting the fact that the row and column associated with 2416a zero diagonal usually have structural ones placed symmetrically. For 2417this routine to work, the structural ones must be exactly equal to either 2418one or negative one. This routine should be used only on modified node ad- 2419mittance matrices and must be executed after the matrix has been built but 2420before spScale(), spNorm(), spMultiply(), spFactor(), spOrderAndFactor() or 2421spDeleteRowAndCol() are executed. It should be executed for the initial 2422factorization only. 2423 2424void spMNA Preorder( Matrix ) 2425 2426o Argument: 2427 2428 Matrix input (char *) 2429 2430 Pointer to the matrix to be preordered. 2431 2432o Compiler options that must be set for this routine to exist: 2433 MODIFIED NODAL 2434 2435 2436 2437 2438 2439 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453 2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465 2466 June 23, 1988 2467 2468 2469 2470 2471 2472 - 37 - 2473 2474 2475 2476 24774.24: spMultiply() 2478 2479Multiplies Matrix by Solution on the right to find RHS. Assumes matrix has 2480not been factored. This routine can be used as a test to see if solutions 2481are correct. 2482 2483void spMultiply( Matrix, RHS, Solution ) 2484void spMultiply( Matrix, RHS, Solution, iRHS, iSolution ) 2485 2486o Arguments: 2487 2488 Matrix input (char *) 2489 Pointer to the matrix. 2490 2491 RHS output (spREAL[]) 2492 RHS is the right hand side vector. This is what is being solved 2493 for. RHS contains only the real portion of the right-hand side 2494 if spSEPARATED COMPLEX VECTORS is set true. 2495 2496 Solution input (spREAL[]) 2497 Solution is the vector being multiplied by the matrix. Solution 2498 contains only the real portion of that vector if 2499 spSEPARATED COMPLEX VECTORS is set true. 2500 2501 iRHS output (spREAL[]) 2502 iRHS is the imaginary portion of the right hand side. This is 2503 what is being solved for. It is only necessary to supply iRHS if 2504 the matrix is complex and spSEPARATED COMPLEX VECTORS is set 2505 true. 2506 2507 iSolution input (spREAL[]) 2508 iSolution is the imaginary portion of the vector being multiplied 2509 by the matrix. It is only necessary to supply iRHS if the matrix 2510 is complex and spSEPARATED COMPLEX VECTORS is set true. 2511 2512o Compiler options that must be set for this routine to exist: 2513 MULTIPLICATION 2514 2515 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2526 2527 2528 2529 2530 2531 2532 June 23, 1988 2533 2534 2535 2536 2537 2538 - 38 - 2539 2540 2541 2542 2543 25444.25: spMultTransposed() 2545 2546Multiplies transposed Matrix by Solution on the right to find RHS. Assumes 2547matrix has not been factored. This routine can be used as a test to see if 2548solutions are correct. 2549 2550void spMultTransposed( Matrix, RHS, Solution ) 2551void spMultTransposed( Matrix, RHS, Solution, iRHS, iSolution ) 2552 2553o Arguments: 2554 2555 Matrix input (char *) 2556 Pointer to the matrix. 2557 2558 RHS output (spREAL[]) 2559 RHS is the right hand side vector. This is what is being solved 2560 for. RHS contains only the real portion of the right-hand side 2561 if spSEPARATED COMPLEX VECTORS is set true. 2562 2563 Solution input (spREAL[]) 2564 Solution is the vector being multiplied by the matrix. Solution 2565 contains only the real portion of that vector if 2566 spSEPARATED COMPLEX VECTORS is set true. 2567 2568 iRHS output (spREAL[]) 2569 iRHS is the imaginary portion of the right hand side. This is 2570 what is being solved for. It is only necessary to supply iRHS if 2571 the matrix is complex and spSEPARATED COMPLEX VECTORS is set 2572 true. 2573 2574 iSolution input (spREAL[]) 2575 iSolution is the imaginary portion of the vector being multiplied 2576 by the matrix. It is only necessary to supply iRHS if the matrix 2577 is complex and spSEPARATED COMPLEX VECTORS is set true. 2578 2579o Compiler options that must be set for this routine to exist: 2580 MULTIPLICATION 2581 TRANSPOSE 2582 2583 2584 2585 2586 2587 2588 2589 2590 2591 2592 2593 2594 2595 2596 2597 2598 2599 June 23, 1988 2600 2601 2602 2603 2604 2605 - 39 - 2606 2607 2608 2609 2610 26114.26: spNorm() 2612 2613Computes and returns the L-infinity norm of an unfactored matrix. This 2614number is used in computing the condition number of the matrix. It is a 2615fatal error to pass this routine a factored matrix. 2616 2617spREAL spNorm( Matrix ) 2618 2619o Returns: 2620 The largest absolute row sum (the L-infinity norm) of the matrix. 2621 2622o Argument: 2623 2624 Matrix input (char *) 2625 Pointer to the matrix. 2626 2627o Compiler options that must be set for this routine to exist: 2628 CONDITION 2629 2630 2631 2632 2633 26344.27: spOrderAndFactor() 2635 2636This routine chooses a pivot order for the matrix and factors it into LU 2637form. It handles both the initial factorization and subsequent factoriza- 2638tions when a reordering or threshold pivoting is desired. This is handled 2639in a manner that is transparent to the user. 2640 2641int spOrderAndFactor( Matrix, RHS, Threshold, AbsoluteThreshold, DiagPivot- 2642ing ) 2643 2644o Returned: 2645 The error code is returned. Possible errors are listed below. 2646 2647o Arguments: 2648 2649 Matrix input (char *) 2650 Pointer to matrix to be factored. 2651 2652 RHS input (spREAL[]) 2653 Representative RHS vector that is used to determine pivoting 2654 order when the right-hand side vector is sparse. If a term in 2655 RHS is zero, it is assumed that it will usually be zero. Con- 2656 versely, a nonzero term in RHS indicates that the term will often 2657 be nonzero. If RHS is a NULL pointer then the right-hand side 2658 vector is assumed to be full and it is not used when determining 2659 the pivoting order. 2660 2661 Threshold input (spREAL) 2662 This is the pivot threshold, which should be between zero and 2663 one. If it is one then the pivoting method becomes complete 2664 2665 2666 2667 June 23, 1988 2668 2669 2670 2671 2672 2673 - 40 - 2674 2675 2676 pivoting, which is very slow and tends to fill up the matrix. If 2677 it is set close to zero the pivoting method becomes strict Mar- 2678 kowitz with no threshold. The pivot threshold is used to elim- 2679 inate pivot candidates that would cause excessive element growth 2680 if they were used. Element growth is the cause of roundoff 2681 error, which can occur even in well-conditioned matrices. Set- 2682 ting the threshold large will reduce element growth and roundoff 2683 error, but setting it too large will cause execution time to be 2684 excessive and will result in a large number of fill-ins. If this 2685 occurs, accuracy can actually be degraded because of the large 2686 number of operations required on the matrix due to the large 2687 number of fill-ins. A good value for diagonal pivoting seems to 2688 be 0.001 while a good value for complete pivoting appears to be 2689 0.1. The default is chosen by giving a value larger than one or 2690 less than or equal to zero. Once the pivot threshold is set, the 2691 value becomes the new default for later calls to spOrderAndFac- 2692 tor. The threshold value should be increased and the matrix re- 2693 solved if growth is found to be excessive. Changing the pivot 2694 threshold does not improve performance on matrices where growth 2695 is low, as is often the case with ill-conditioned matrices. The 2696 default value of Threshold was choosen for use with nearly diago- 2697 nally dominant matrices such as node- and modified-node admit- 2698 tance matrices. For these matrices it is usually best to use 2699 diagonal pivoting. For matrices without a strong diagonal, it is 2700 usually best to use a larger threshold, such as 0.01 or 0.1. 2701 2702 AbsoluteThreshold input (spREAL) 2703 The absolute magnitude an element must have to be considered as a 2704 pivot candidate, except as a last resort. This number should be 2705 set significantly smaller than the smallest diagonal element that 2706 is is expected to be placed in the matrix. If there is no rea- 2707 sonable prediction for the lower bound on these elements, then 2708 AbsoluteThreshold should be set to zero. AbsoluteThreshold is 2709 used to reduce the possibility of choosing as a pivot an element 2710 that has suffered heavy cancellation and as a result mainly con- 2711 sists of roundoff error. Note that if AbsoluteThreshold is set 2712 too large, it could drastically increase the time required to 2713 factor and solve the matrix. AbsoluteThreshold should be nonne- 2714 gative. If no element in the matrix is larger than Absolu- 2715 teThreshold, the warning spILL CONDITIONED is returned. 2716 2717 DiagPivoting input (int) 2718 A flag indicating that pivot selection should be confined to the 2719 diagonal if possible. If DiagPivoting is nonzero and if 2720 DIAGONAL PIVOTING is enabled pivots will be chosen only from the 2721 diagonal unless there are no diagonal elements that satisfy the 2722 threshold criteria. Otherwise, the entire reduced submatrix is 2723 searched when looking for a pivot. The diagonal pivoting in 2724 Sparse is efficient and well refined, while the complete pivoting 2725 is not. For symmetric and near symmetric matrices, it is best to 2726 use diagonal pivoting because it results in the best performance 2727 when reordering the matrix and when factoring the matrix without 2728 ordering. If there is a considerable amount of nonsymmetry in 2729 the matrix, then complete pivoting may result in a better 2730 2731 2732 2733 June 23, 1988 2734 2735 2736 2737 2738 2739 - 41 - 2740 2741 2742 equation ordering simply because there are more pivot candidates 2743 to choose from. A better ordering results in faster subsequent 2744 factorizations. However, the initial pivot selection process 2745 takes considerably longer for complete pivoting. 2746 2747o Possible errors: 2748 spNO MEMORY 2749 spSINGULAR 2750 spILL CONDITIONED 2751 2752 2753 2754 2755 27564.28: spPartition() 2757 2758This routine determines the cost to factor each row using both direct and 2759indirect addressing and decides, on a row-by-row basis, which addressing 2760mode is fastest. This information is used in spFactor() to speed the fac- 2761torization. 2762 2763When factoring a previously ordered matrix using spFactor(), fISparse 2764operates on a row-at-a-time basis. For speed, on each step, the row being 2765updated is copied into a full vector and the operations are performed on 2766that vector. This can be done one of two ways, either using direct ad- 2767dressing or indirect addressing. Direct addressing is fastest when the ma- 2768trix is relatively dense and indirect addressing is best when the matrix is 2769quite sparse. The user selects the type of partition used with Mode. If 2770Mode is set to spDIRECT PARTITION, then the all rows are placed in the 2771direct addressing partition. Similarly, if Mode is set to 2772spINDIRECT PARTITION, then the all rows are placed in the indirect address- 2773ing partition. By setting Mode to spAUTO PARTITION, the user allows Sparse 2774to select the partition for each row individually. spFactor() generally 2775runs faster if Sparse is allowed to choose its own partitioning, however 2776choosing a partition is expensive. The time required to choose a partition 2777is of the same order of the cost to factor the matrix. If you plan to fac- 2778tor a large number of matrices with the same structure, it is best to let 2779Sparse choose the partition. Otherwise, you should choose the partition 2780based on the predicted density of the matrix. By default (i.e., if spPar- 2781tition() is never called), Sparse chooses the partition for each row indi- 2782vidually. 2783 2784void spPartition( Matrix, Mode ) 2785 2786o Arguments: 2787 2788 Matrix input (char *) 2789 Pointer to matrix to be partitioned. 2790 2791 Mode input (int) 2792 Mode must be one of three special codes: spDIRECT PARTITION, 2793 spINDIRECT PARTITION, or spAUTO PARTITION. 2794 2795 2796 2797 2798 2799 2800 June 23, 1988 2801 2802 2803 2804 2805 2806 - 42 - 2807 2808 2809 2810 2811 28124.29: spPrint() 2813 2814Formats and send the matrix to standard output. Some elementary statistics 2815are also output. The matrix is output in a format that is readable by peo- 2816ple. This routine should not be used on large matrices. 2817 2818void spPrint( Matrix, PrintReordered, Data, Header ) 2819 2820o Arguments: 2821 2822 Matrix input (char *) 2823 Pointer to matrix to be printed. 2824 2825 PrintReordered input (int) 2826 Indicates whether the matrix should be printed out in its origi- 2827 nal form, as input by the user, or whether it should be printed 2828 in its reordered form, as used internally by the matrix routines. 2829 A zero indicates that the matrix should be printed as inputed, a 2830 one indicates that it should be printed reordered. 2831 2832 Data input (int) 2833 Boolean flag that when false indicates that output should be 2834 compressed such that only the existence of an element should be 2835 indicated rather than giving the actual value. Thus 10 times as 2836 many elements can be printed on a row. A zero indicates that the 2837 matrix should be printed compressed. A one signifies that the 2838 matrix should be printed in all its glory. 2839 2840 Header input (int) 2841 A flag indicating that extra information should be printed, such 2842 as row and column numbers. 2843 2844o Compiler options that must be set for this routine to exist: 2845 DOCUMENTATION 2846 2847 2848 2849 2850 2851 2852 2853 2854 2855 2856 2857 2858 2859 2860 2861 2862 2863 2864 2865 2866 2867 June 23, 1988 2868 2869 2870 2871 2872 2873 - 43 - 2874 2875 2876 2877 2878 28794.30: spPseudoCondition() 2880 2881Computes the magnitude of the ratio of the largest to the smallest pivots. 2882This quantity is an indicator of ill-conditioning in the matrix. If this 2883ratio is large, and if the matrix is scaled such that uncertainties in the 2884right-hand side vector and the matrix entries are equilibrated, then the 2885matrix is ill-conditioned. However, a small ratio does not necessarily im- 2886ply that the matrix is well-conditioned. This routine must only be used 2887after a matrix has been factored by spOrderAndFactor() or spFactor() and 2888before it is cleared by spClear() or spInitialize(). The pseudocondition 2889is faster to compute than the condition number calculated by spCondition(), 2890but is not as informative. 2891 2892spREAL spPseudoCondition( Matrix ) 2893 2894o Returns: 2895 The magnitude of the ratio of the largest to smallest pivot used dur- 2896 ing previous factorization. If the matrix was singular, zero is re- 2897 turned. 2898 2899o Argument: 2900 2901 Matrix input (char *) 2902 Pointer to matrix. 2903 2904o Compiler options that must be set for this routine to exist: 2905 PSEUDOCONDITION 2906 2907 2908 2909 2910 2911 2912 2913 2914 2915 2916 2917 2918 2919 2920 2921 2922 2923 2924 2925 2926 2927 2928 2929 2930 2931 2932 2933 2934 June 23, 1988 2935 2936 2937 2938 2939 2940 - 44 - 2941 2942 2943 2944 2945 29464.31: spRoundoff() 2947 2948Returns a bound on the magnitude of the largest element in E = A-LU, where 2949E represents error in the matrix resulting from roundoff during the factor- 2950ization. 2951 2952spREAL spRoundoff( Matrix, Rho ) 2953 2954o Returns: 2955 Returns a bound on the magnitude of the largest element in E = A-LU. 2956 2957o Arguments: 2958 2959 Matrix input (char *) 2960 Pointer to matrix. Matrix must be factored. 2961 2962 Rho input (spREAL) 2963 The bound on the magnitude of the largest element in any of the 2964 reduced submatrices. This is the number computed by the function 2965 spLargestElement() when given a factored matrix. If this number 2966 is negative, the bound will be computed automatically. 2967 2968o Compiler options that must be set for this routine to exist: 2969 STABILITY 2970 2971 2972 2973 2974 2975 2976 2977 2978 2979 2980 2981 2982 2983 2984 2985 2986 2987 2988 2989 2990 2991 2992 2993 2994 2995 2996 2997 2998 2999 3000 3001 June 23, 1988 3002 3003 3004 3005 3006 3007 - 45 - 3008 3009 3010 3011 30124.32: spScale() 3013 3014This function scales the matrix to enhance the possibility of finding a 3015good pivoting order. Note that scaling enhances accuracy of the solution 3016only if it affects the pivoting order, so it only makes sense to scale the 3017matrix before spOrderAndFactor(). There are several things to take into 3018account when choosing the scale factors. First, the scale factors are 3019directly multiplied times the elements in the matrix. To prevent roundoff, 3020each scale factor should be equal to an integer power of the number base of 3021the machine. Since most machines operate in base two, scale factors should 3022be a power of two. Second, the matrix should be scaled such that the ma- 3023trix of element uncertainties is equilibrated. Third, this function multi- 3024plies the scale factors times the elements, so if one row tends to have un- 3025certainties 1000 times smaller than the other rows, then its scale factor 3026should be 1024, not 1/1024. Fourth, to save time, this function does not 3027scale rows or columns if their scale factors are equal to one. Thus, the 3028scale factors should be normalized to the most common scale factor. Rows 3029and columns should be normalized separately. For example, if the size of 3030the matrix is 100 and 10 rows tend to have uncertainties near 1e-6 and the 3031remaining 90 have uncertainties near 1e-12, then the scale factor for the 303210 should be 1/1,048,576 and the scale factors for the remaining 90 should 3033be 1. Fifth, since this routine directly operates on the matrix, it is 3034necessary to apply the scale factors to the RHS and Solution vectors. It 3035may be easier to simply use spOrderAndFactor() on a scaled matrix to choose 3036the pivoting order, and then throw away the matrix. Subsequent factoriza- 3037tions, performed with spFactor(), will not need to have the RHS and Solu- 3038tion vectors descaled. 3039 3040void spScale( Matrix, RHS ScaleFactors, SolutionScaleFactors ) 3041 3042o Arguments: 3043 3044 Matrix input (char *) 3045 Pointer to the matrix to be scaled. 3046 3047 RHS ScaleFactors input (spREAL[]) 3048 The array of RHS scale factors. These factors scale the rows. 3049 All scale factors are real-valued. 3050 3051 SolutionScaleFactors input (spREAL[]) 3052 The array of Solution scale factors. These factors scale the 3053 columns. All scale factors are real-valued. 3054 3055o Compiler options that must be set for this routine to exist: 3056 SCALING 3057 3058 3059 3060 3061 3062 3063 3064 3065 3066 3067 June 23, 1988 3068 3069 3070 3071 3072 3073 - 46 - 3074 3075 3076 3077 3078 30794.33: spSetComplex() 3080 3081The type of the matrix may then be toggled back and forth between complex 3082and real. This function changes the type of matrix to complex. For the 3083matrix to be set complex, the compiler option spCOMPLEX must be set true. 3084 3085void spSetComplex( Matrix ) 3086 3087o Argument: 3088 3089 Matrix input (char *) 3090 3091 The matrix that is to be to be complex. 3092 3093 3094 3095 3096 30974.34: spSetReal() 3098 3099The type of the matrix may then be toggled back and forth between complex 3100and real. This function changes the type of matrix to real. For the ma- 3101trix to be set real, the compiler option REAL must be set true. 3102 3103void spSetReal( Matrix ) 3104 3105o Argument: 3106 3107 Matrix input (char *) 3108 The matrix that is to be real. 3109 3110 3111 3112 3113 3114 3115 3116 3117 3118 3119 3120 3121 3122 3123 3124 3125 3126 3127 3128 3129 3130 3131 3132 3133 3134 3135 June 23, 1988 3136 3137 3138 3139 3140 3141 - 47 - 3142 3143 3144 3145 3146 31474.35: spSolve() 3148 3149Performs the forward and backward elimination to find the unknown Solution 3150vector from RHS and the factored matrix. 3151 3152void spSolve( Matrix, RHS, Solution ) 3153void spSolve( Matrix, RHS, Solution, iRHS, iSolution ) 3154 3155o Arguments: 3156 3157 Matrix input (char *) 3158 Pointer to matrix. 3159 3160 RHS input (spREAL[]) 3161 RHS is the input data array, the right-hand side vector. RHS con- 3162 tains only the real portion of the right-hand side vector if 3163 spSEPARATED COMPLEX VECTORS is set true. RHS is undisturbed and 3164 may be reused for other solves. 3165 3166 Solution output (spREAL[]) 3167 Solution is the output data array, the unknown vector. This rou- 3168 tine is constructed such that RHS and Solution can be the same 3169 array. Solution contains only the real portion of the unknown 3170 vector if spSEPARATED COMPLEX VECTORS is set true. 3171 3172 iRHS input (spREAL[]) 3173 iRHS is the imaginary portion of the input data array, the 3174 right-hand side vector. This data is undisturbed and may be 3175 reused for other solves. This argument is unnecessary if the ma- 3176 trix is real or spSEPARATED COMPLEX VECTORS is set false. 3177 3178 iSolution output (spREAL[]) 3179 iSolution is the imaginary portion of the output data array. 3180 This routine is constructed such that iRHS and iSolution can be 3181 the same array. This argument is unnecessary if the matrix is 3182 real or spSEPARATED COMPLEX VECTORS is set false. 3183 3184 3185 3186 3187 3188 3189 3190 3191 3192 3193 3194 3195 3196 3197 3198 3199 3200 3201 3202 June 23, 1988 3203 3204 3205 3206 3207 3208 - 48 - 3209 3210 3211 3212 3213 32144.36: spSolveTransposed() 3215 3216Performs the forward and backward elimination to find the unknown Solution 3217vector from RHS and the transposed factored matrix. This routine is useful 3218when performing sensitivity analysis on a circuit using the adjoint method. 3219 3220void spSolveTransposed( Matrix, RHS, Solution ) 3221void spSolveTransposed( Matrix, RHS, Solution, iRHS, iSolution ) 3222 3223o Arguments: 3224 3225 Matrix input (char *) 3226 Pointer to matrix. 3227 3228 RHS input (spREAL[]) 3229 RHS is the input data array, the right-hand side vector. RHS 3230 contains only the real portion of the right-hand side vector if 3231 spSEPARATED COMPLEX VECTORS is set true. RHS is undisturbed and 3232 may be reused for other solves. 3233 3234 Solution output (spREAL[]) 3235 Solution is the output data array, the unknown vector. This rou- 3236 tine is constructed such that RHS and Solution can be the same 3237 array. Solution contains only the real portion of the unknown 3238 vector if spSEPARATED COMPLEX VECTORS is set true. 3239 3240 iRHS input (spREAL[]) 3241 iRHS is the imaginary portion of the input data array, the 3242 right-hand side vector. This data is undisturbed and may be 3243 reused for other solves. This parameter is unnecessary if the 3244 matrix is real or spSEPARATED COMPLEX VECTORS is set false. 3245 3246 iSolution output (spREAL[]) 3247 iSolution is the imaginary portion of the output data array. 3248 This routine is constructed such that iRHS and iSolution can be 3249 the same array. This parameter is unnecessary if the matrix is 3250 real or spSEPARATED COMPLEX VECTORS is set false. 3251 3252o Compiler options that must be set for this routine to exist: 3253 TRANSPOSE 3254 3255 3256 3257 3258 3259 3260 3261 3262 3263 3264 3265 3266 3267 3268 3269 June 23, 1988 3270 3271 3272 3273 3274 3275 - 49 - 3276 3277 3278 3279 3280 32814.37: spStripFills() 3282 3283spStripFills() strips all accumulated fill-ins from a matrix. This is 3284often a useful thing to do before reordering a matrix to help insure that 3285subsequent factorizations will be as efficient as possible. 3286 3287void spStripFills( Matrix ) 3288 3289o Argument: 3290 3291 Matrix input (char *) 3292 The matrix to be stripped. 3293 3294o Compiler options that must be set for this routine to exist: 3295 STRIP 3296 3297 3298 3299 3300 33014.38: spWhereSingular() 3302 3303This function returns the row and column number where the matrix was 3304detected as singular or where a zero pivot was found. 3305 3306void spWhereSingular( Matrix, Row, Col ) 3307 3308o Arguments: 3309 3310 Matrix input (char *) 3311 Pointer to matrix. 3312 3313 Row output (int *) 3314 The row number. 3315 3316 Row output (int *) 3317 The column number. 3318 3319 3320 3321 3322 3323 3324 3325 3326 3327 3328 3329 3330 3331 3332 3333 3334 3335 3336 3337 June 23, 1988 3338 3339 3340 3341 3342 3343 - 50 - 3344 3345 33465: MACRO FUNCTIONS 3347These macro functions are used to quickly enter data into the matrix using 3348pointers. These pointers are originally acquired by the user from 3349spGetElement(), spGetAdmittance(), spGetQuad(), and spGetOnes() during the 3350initial loading of the matrix. These macros work correctly even if the 3351elements they are to add data to are in row or column zero. 3352 3353 The macros reside in the file spExports.h. To use them, this file 3354must be included in the file of the calling routine and that routine must 3355be written in C. 3356 3357 33585.1: spADD REAL ELEMENT() 3359 3360Macro function that adds a real value to an element in the matrix by a 3361pointer. 3362 3363spADD REAL ELEMENT( pElement , Real ) 3364 3365o Arguments: 3366 3367 pElement input (spREAL *) 3368 A pointer to the element to which Real is to be added. 3369 3370 Real input (spREAL) 3371 The real value that is to be added to the element. 3372 3373 3374 3375 3376 33775.2: spADD IMAG ELEMENT() 3378 3379Macro function that adds a imaginary value to an element in the matrix by a 3380pointer. 3381 3382spADD IMAG ELEMENT( pElement , Imag ) 3383 3384o Arguments: 3385 3386 pElement input (spREAL *) 3387 A pointer to the element to which Imag is to be added. 3388 3389 Imag input (spREAL) 3390 The imaginary value that is to be added to the element. 3391 3392 3393 3394 3395 3396 3397 3398 3399 3400 3401 3402 3403 3404 June 23, 1988 3405 3406 3407 3408 3409 3410 - 51 - 3411 3412 3413 3414 3415 34165.3: spADD COMPLEX ELEMENT() 3417 3418Macro function that adds a complex value to an element in the matrix by a 3419pointer. 3420 3421spADD COMPLEX ELEMENT( pElement, Real, Imag ) 3422 3423o Arguments: 3424 3425 pElement input (spREAL *) 3426 A pointer to the element to which Real and Imag are to be added. 3427 3428 Real input (spREAL) 3429 The real value that is to be added to the element. 3430 3431 Imag input (spREAL) 3432 The imaginary value that is to be added to the element. 3433 3434 3435 3436 3437 34385.4: spADD REAL QUAD() 3439 3440Macro that adds a real value to the four elements specified by Template. 3441The value is added to the first two elements in Template, and subtracted 3442from the last two. 3443 3444spADD REAL QUAD( Template, Real ) 3445 3446o Arguments: 3447 3448 Template input (struct spTemplate) 3449 Data structure containing the pointers to four matrix elements. 3450 3451 Real input (spREAL) 3452 Real value to be added to the elements. 3453 3454 3455 3456 3457 3458 3459 3460 3461 3462 3463 3464 3465 3466 3467 3468 3469 3470 3471 3472 June 23, 1988 3473 3474 3475 3476 3477 3478 - 52 - 3479 3480 3481 3482 3483 34845.5: spADD IMAG QUAD() 3485 3486Macro that adds an imaginary value to the four elements specified by Tem- 3487plate. The value is added to the first two elements in Template, and sub- 3488tracted from the last two. 3489 3490spADD IMAG QUAD( Template, Imag ) 3491 3492o Arguments: 3493 3494 Template input (struct spTemplate) 3495 Data structure containing the pointers to four matrix elements. 3496 3497 Imag input (spREAL) 3498 Imaginary value to be added to the elements. 3499 3500 3501 3502 3503 35045.6: spADD COMPLEX QUAD() 3505 3506Macro that adds a complex value to the four elements specified by Template. 3507The value is added to the first two elements in Template, and subtracted 3508from the last two. 3509 3510spADD COMPLEX QUAD( Template, Real, Imag ) 3511 3512o Arguments: 3513 3514 Template input (struct spTemplate) 3515 Data structure containing the pointers to four matrix elements. 3516 3517 Real input (spREAL) 3518 Real value to be added to the elements. 3519 3520 Imag input (spREAL) 3521 Imaginary value to be added to the elements. 3522 3523 3524 3525 3526 3527 3528 3529 3530 3531 3532 3533 3534 3535 3536 3537 3538 3539 3540 June 23, 1988 3541 3542 3543 3544 3545 3546 - 53 - 3547 3548 35496: CONFIGURING SPARSE 3550 3551 Sparse has a extensive set of options and parameters that can be set 3552at compile time to alter the personality of the program. They also are 3553used to eliminate routines that are not needed so as to reduce the amount 3554of memory required to hold the object code. These options and parameters 3555consist of macros definitions and are contained in the file spConfig.h. To 3556configure Sparse, spConfig.h must be edited and then Sparse must be recom- 3557piled. 3558 3559 Some terminology should be defined. The Markowitz row count is the 3560number of non-zero elements in a row excluding the one being considered as 3561pivot. There is one Markowitz row count for every row. The Markowitz 3562column count is defined similarly for columns. The Markowitz product for 3563an element is the product of its row and column counts. It is a measure of 3564how much work would be required on the next step of the factorization if 3565that element were chosen to be pivot. A small Markowitz product is desir- 3566able. For a more detailed explanation, see Kundert [kundert86]. 3567 3568 35696.1: Sparse Options 3570 3571REAL 3572 3573This specifies that the routines are expected to handle real systems of 3574equations. The routines can be compiled to handle both real and complex 3575systems at the same time, but there is a slight speed and memory advantage 3576if the routines are complied to handle only real systems of equations. 3577 3578 3579spCOMPLEX 3580 3581This specifies that the routines will be complied to handle complex systems 3582of equations. 3583 3584 3585EXPANDABLE 3586 3587Setting this compiler flag true makes the matrix expandable before it has 3588been reordered. If the matrix is expandable, then if an element is added 3589that would be considered out of bounds in the current matrix, the size of 3590the matrix is increased to hold that element. As a result, the size of the 3591matrix need not be known before the matrix is built. The matrix can be al- 3592located with size zero and expanded. It is possible to expand the size of 3593a matrix after it is been reordered if TRANSLATE and EXPANDABLE are both 3594set true. 3595 3596 3597 3598 3599 3600 3601 3602 3603 3604 3605 3606 3607 June 23, 1988 3608 3609 3610 3611 3612 3613 - 54 - 3614 3615 3616 3617TRANSLATE 3618 3619This option allows the set of external row and column numbers to be non- 3620packed. In other words, the row and column numbers need not be contiguous. 3621The priced paid for this flexibility is that when TRANSLATE is set true, 3622the time required to initially build the matrix will be greater because the 3623external row and column number must be translated into internal 3624equivalents. This translation brings about other benefits though. First, 3625the spGetElement(), spGetAdmittance(), spGetQuad(), and spGetOnes() rou- 3626tines may be used after the matrix has been factored. Further, elements, 3627and even rows and columns, may be added to the matrix, and rows and columns 3628may be deleted from the matrix, after it has been reordered. Note that 3629when the set of row and column number is not a packed set, neither are the 3630RHS and Solution vectors. Thus the size of these vectors must be at least 3631as large as the external size, which is the value of the largest given row 3632or column numbers. 3633 3634 3635INITIALIZE 3636 3637Causes the spInitialize(), spGetInitInfo(), and spInstallInitInfo() rou- 3638tines to be compiled. These routines allow the user to store and read one 3639pointer in each nonzero element in the matrix. spInitialize() then calls a 3640user specified function for each structural nonzero in the matrix, and in- 3641cludes this pointer as well as the external row and column numbers as argu- 3642ments. This allows the user to write custom matrix and right-hand side 3643vector initialization routines. 3644 3645 3646DIAGONAL PIVOTING 3647 3648Many matrices, and in particular node- and modified-node admittance ma- 3649trices, tend to be nearly symmetric and nearly diagonally dominant. For 3650these matrices, it is a good idea to select pivots from the diagonal. With 3651this option enabled, this is exactly what happens, though if no satisfacto- 3652ry pivot can be found on the diagonal, an off-diagonal pivot will be used. 3653If this option is disabled, Sparse does not preferentially search the diag- 3654onal. Because of this, Sparse has a wider variety of pivot candidates 3655available, and so presumably fewer fill-ins will be created. However, the 3656initial pivot selection process will take considerably longer. If working 3657with node admittance matrices, or other matrices with a strong diagonal, it 3658is probably best to use DIAGONAL PIVOTING for two reasons. First, accuracy 3659will be better because pivots will be chosen from the large diagonal ele- 3660ments, thus reducing the chance of growth and hence, roundoff. Second, a 3661near optimal ordering will be chosen quickly. If the class of matrices you 3662are working with does not have a strong diagonal, do not use 3663DIAGONAL PIVOTING, but consider using a larger threshold. When 3664DIAGONAL PIVOTING is turned off, the following options and constants are 3665not used: MODIFIED MARKOWITZ, MAX MARKOWITZ TIES, and TIES MULTIPLIER. 3666 3667 3668 3669 3670 3671 3672 3673 3674 June 23, 1988 3675 3676 3677 3678 3679 3680 - 55 - 3681 3682 3683 3684ARRAY OFFSET 3685 3686This determines whether arrays start at an index of zero or one. This op- 3687tion is necessitated by the fact that standard C convention dictates that 3688arrays begin with an index of zero but the standard mathematic convention 3689states that arrays begin with an index of one. So if you prefer to start 3690your arrays with zero, or you're calling Sparse from some other programming 3691language, use an ARRAY OFFSET of 0. Otherwise, use an ARRAY OFFSET of 1. 3692Note that if you use an offset of one, the arrays that you pass to Sparse 3693must have an allocated length of one plus the external size of the matrix. 3694ARRAY OFFSET must be either 0 or 1, no other offsets are valid. 3695 3696 3697spSEPARATED COMPLEX VECTORS 3698 3699This specifies the format for complex vectors. If this is set false then a 3700complex vector is made up of one double sized array of spREALs in which the 3701real and imaginary numbers are placed alternately in the array. In other 3702words, the first entry would be Complex[1].Real, then comes 3703Complex[1].Imag, then Complex[2].Real, etc. If spSEPARATED COMPLEX VECTORS 3704is set true, then each complex vector is represented by two arrays of 3705spREALs, one with the real terms, the other with the imaginary. 3706 3707 3708MODIFIED MARKOWITZ 3709 3710This specifies that the modified Markowitz method of pivot selection is to 3711be used. The modified Markowitz method differs from standard Markowitz in 3712two ways. First, under modified Markowitz, the search for a pivot can be 3713terminated early if a adequate (in terms of sparsity) pivot candidate is 3714found. Thus, when using modified Markowitz, the initial factorization can 3715be faster, but at the expense of a suboptimal pivoting order that may slow 3716subsequent factorizations. The second difference is in the way modified 3717Markowitz breaks Markowitz ties. When two or more elements are pivot can- 3718didates and they all have the same Markowitz product, then the tie is bro- 3719ken by choosing the element that is best numerically. The numerically best 3720element is the one with the largest ratio of its magnitude to the magnitude 3721of the largest element in the same column, excluding itself. The modified 3722Markowitz method results in marginally better accuracy. 3723 3724 3725DELETE 3726 3727This specifies that the spDeleteRowAndCol() routine should be compiled. 3728Note that for this routine to be compiled, both DELETE and TRANSLATE should 3729be set true. 3730 3731 3732STRIP 3733 3734This specifies that the spStripFills() routine should be compiled. 3735 3736 3737 3738 3739 3740 3741 3742 June 23, 1988 3743 3744 3745 3746 3747 3748 - 56 - 3749 3750 3751 3752MODIFIED NODAL 3753 3754This specifies that the spMNA Preorder(), the routine that preorders modi- 3755fied node admittance matrices, should be compiled. This routine results in 3756greater speed and accuracy if used with this type of matrix. 3757 3758 3759QUAD ELEMENT 3760 3761This specifies that the routines that allow four related elements to be en- 3762tered into the matrix at once should be compiled. The routines affected by 3763QUAD ELEMENT are spGetAdmittance(), spGetQuad(), and spGetOnes(). 3764 3765 3766TRANSPOSE 3767 3768This specifies that spSolveTranspose() and perhaps spMultTransposed(), 3769which operate on the matrix as if it was transposed, should be compiled. 3770 3771SCALING 3772 3773This specifies that the routine that performs scaling on the matrix should 3774be complied. Scaling is not strongly supported. The routine to scale the 3775matrix is provided, but no routines are provided to scale and descale the 3776RHS and Solution vectors. It is suggested that if scaling is desired, it 3777only be performed when the pivot order is being chosen, which is done in 3778spOrderAndFactor(). This, and when the condition number of the matrix is 3779calculated with spCondition(), are the only times scaling has an effect. 3780The scaling may then either be removed from the solution by the user or the 3781scaled factored matrix may simply be thrown away. 3782 3783 3784DOCUMENTATION 3785 3786This specifies that routines that are used to document the matrix, 3787spPrint(), spFileMatrix(), spFileVector(), and spFileStats(), should be 3788compiled. 3789 3790 3791DETERMINANT 3792 3793This specifies that the spDeterminant() routine should be complied. 3794 3795 3796STABILITY 3797 3798This specifies that spLargestElement() and spRoundoff() should be compiled. 3799These routines are used to check the stability (and hence the quality of 3800the pivoting) of the factorization by computing a bound on the size of the 3801element is the matrix E = A-LU. If this bound is very high after applying 3802spOrderAndFactor(), then the pivot threshold should be raised. If the 3803bound increases greatly after using spFactor(), then the matrix should 3804probably be reordered. 3805 3806 3807 3808 3809 3810 June 23, 1988 3811 3812 3813 3814 3815 3816 - 57 - 3817 3818 3819 3820CONDITION 3821 3822This specifies that spCondition() and spNorm(), the code that computes a 3823good estimate of the condition number of the matrix, should be compiled. 3824 3825 3826PSEUDOCONDITION 3827 3828This specifies that spPseudoCondition(), the code that computes a crude and 3829easily fooled indicator of the ill-conditioning in the matrix, should be 3830compiled. 3831 3832 3833MULTIPLICATION 3834 3835This specifies that spMultiply() and perhaps spMultTransposed(), the rou- 3836tines that multiply an unfactored matrix by a vector, should be compiled. 3837 3838 3839FORTRAN 3840 3841This specifies that the FORTRAN interface to Sparse1.3 should be compiled. 3842The ARRAY OFFSET option should be set to NO when interfacing to FORTRAN 3843programs. 3844 3845 3846DEBUG 3847 3848This specifies that additional error checking should be compiled. The type 3849of errors checked are those that are common when the matrix routines are 3850first integrated into a user's program. Once the routines have been in- 3851tegrated in and are running smoothly, this option should be turned off. 3852With DEBUG enabled, Sparse is very defensive. If a Sparse routine is 3853called improperly, a message will be printed describing the file and line 3854number where the error was found and execution is aborted. One thing that 3855Sparse is particularly picky about is calling certain functions after an 3856error has occurred. If an error has occurred, do not call 3857spMNA Preorder(), spScale(), spOrderAndFactor(), spFactor(), spSolve(), or 3858spSolveTransposed() until the error has been cleared by spClear() or spIni- 3859tialize(). 3860 3861 3862spCOMPATIBILITY 3863 3864This specifies that Sparse1.3 should be configured to be upward compatible 3865from Sparse1.2. This option is not suggested for use in new software. 3866Sparse1.3, when configured to be compatible with Sparse1.2, is not com- 3867pletely compatible. In particular, if recompiling the calling program, it 3868is necessary to change the names of the Sparse1.2 include files. This op- 3869tion will go away on any future version of Sparse. 3870 3871 3872 3873 3874 3875 3876 3877 3878 3879 June 23, 1988 3880 3881 3882 3883 3884 3885 - 58 - 3886 3887 38886.2: Sparse Constants 3889 3890 These constants are used throughout the sparse matrix routines. They 3891should be set to suit the type of matrices being solved. 3892 3893 3894DEFAULT THRESHOLD 3895 3896The threshold used if the user enters an invalid threshold. Also the 3897threshold used by spFactor() when calling spOrderAndFactor(). The default 3898threshold should not be less than or equal to zero nor larger than one. 3899 3900 3901DIAG PIVOTING AS DEFAULT 3902 3903This indicates whether spOrderAndFactor() should use diagonal pivoting as 3904default. This issue only arises when spOrderAndFactor() is called from 3905spFactor(). 3906 3907 3908SPACE FOR ELEMENTS 3909 3910This number multiplied by the size of the matrix equals the number of ele- 3911ments for which memory is initially allocated in spCreate(). 3912 3913 3914SPACE FOR FILL INS 3915 3916This number multiplied by the size of the matrix equals the number of ele- 3917ments for which memory is initially allocated and specifically reserved for 3918fill-ins in spCreate(). 3919 3920 3921ELEMENTS PER ALLOCATION 3922 3923The number of matrix elements requested from the malloc utility on each 3924call to it. Setting this value greater than one reduces the amount of 3925overhead spent in this system call. 3926 3927 3928MINIMUM ALLOCATED SIZE 3929 3930The minimum allocated size of a matrix. Note that this does not limit the 3931minimum size of a matrix. This just prevents having to resize a matrix 3932many times if the matrix is expandable, large and allocated with an es- 3933timated size of zero. This number must not be less than one. 3934 3935 3936EXPANSION FACTOR 3937 3938The minimum increase in the allocated size of the matrix when it is expand- 3939ed. This number must be greater than one but shouldn't be much larger than 3940two. 3941 3942 3943 3944 3945 3946 3947 3948 3949 June 23, 1988 3950 3951 3952 3953 3954 3955 - 59 - 3956 3957 3958 3959MAX MARKOWITZ TIES 3960 3961This number is used for two slightly different things, both of which relate 3962to the search for the best pivot. First, it is the maximum number of ele- 3963ments that are Markowitz tied that will be sifted through when trying to 3964find the one that is numerically the best. Second, it creates an upper 3965bound on how large a Markowitz product can be before it eliminates the pos- 3966sibility of early termination of the pivot search. In other words, if the 3967product of the smallest Markowitz product yet found and TIES MULTIPLIER is 3968greater than MAX MARKOWITZ TIES, then no early termination takes place. 3969Set MAX MARKOWITZ TIES to some small value if no early termination of the 3970pivot search is desired. An array of spREALs is allocated of size 3971MAX MARKOWITZ TIES so it must be positive and shouldn't be too large. 3972 3973 3974TIES MULTIPLIER 3975 3976Specifies the number of Markowitz ties that are allowed to occur before the 3977search for the pivot is terminated early. Set to some large value if no 3978early termination of the pivot search is desired. This number is multi- 3979plied by the Markowitz product to determine how many ties are required for 3980early termination. This means that more elements will be searched before 3981early termination if a large number of fill-ins could be created by accept- 3982ing what is currently considered the best choice for the pivot. Setting 3983this number to zero effectively eliminates all pivoting, which should be 3984avoided. This number must be positive. 3985 3986 3987DEFAULT PARTITION 3988 3989Which partition mode is used by spPartition() as default. Possibilities 3990include: 3991 3992 spDIRECT PARTITION - each row used direct addressing, best for a few 3993 relatively dense matrices. 3994 3995 spINDIRECT PARTITION - each row used indirect addressing, best for a 3996 few very sparse matrices. 3997 3998 spAUTO PARTITION - direct or indirect addressing is chosen on a row- 3999 by-row basis, carries a large overhead, but speeds up both dense 4000 and sparse matrices, best if there is a large number of matrices 4001 that can use the same ordering. 4002 4003 4004PRINTER WIDTH 4005 4006Gives the number of characters printable in one page width. Set to 80 for 4007terminals and 132 for line printers. 4008 4009 4010 4011 4012 4013 4014 4015 4016 4017 June 23, 1988 4018 4019 4020 4021 4022 4023 - 60 - 4024 4025 40266.3: Machine Constants 4027 4028These numbers must be updated when the program is ported to a new machine. 4029 4030 4031MACHINE RESOLUTION 4032 4033This is the smallest positive real double precision number e such that 40341 + e = 1. 4035 4036 4037LARGEST REAL 4038 4039The largest positive real number representable by a double. 4040 4041 4042SMALLEST REAL 4043 4044The smallest positive real number representable by a double. 4045 4046 4047LARGEST SHORT INTEGER 4048 4049The largest positive integer representable by a short. 4050 4051 4052LARGEST LONG INTEGER 4053 4054The largest positive integer representable by a long. 4055 4056 4057 4058 4059 4060 4061 4062 4063 4064 4065 4066 4067 4068 4069 4070 4071 4072 4073 4074 4075 4076 4077 4078 4079 4080 4081 4082 4083 4084 4085 4086 June 23, 1988 4087 4088 4089 4090 4091 4092 - 61 - 4093 4094 40957: EXPORTS 4096 40977.1: Error Codes 4098 4099 Errors are indicated with a integer error code. Macros definitions 4100for these error codes are set up and placed in the file spMatrix.h. They 4101may be imported into the users program to give readable names to the possi- 4102ble matrix errors. The possible error codes and there corresponding macros 4103are: 4104 4105 4106 4107spOKAY - 0 4108 4109No error has occurred. 4110 4111spSMALL PIVOT - 1 4112 4113When reordering the matrix, no element was found which satisfies the abso- 4114lute threshold criteria. The largest element in the matrix was chosen as 4115pivot. Nonfatal. 4116 4117spZERO DIAG - 2 4118 4119Fatal error. A zero was encountered on the diagonal of the matrix. This 4120does not necessarily imply that the matrix is singular. When this error 4121occurs, the matrix should be reconstructed and factored using spOr- 4122derAndFactor(). 4123 4124spSINGULAR - 3 4125 4126Fatal error. Matrix is singular, so no unique solution exists. 4127 4128spNO MEMORY - 4 4129 4130Fatal error. Indicates that not enough memory is available from the system 4131to handle the matrix. 4132 4133spPANIC - 5 4134 4135Fatal error indicating that the routines are being asked to do something 4136nonsensical or something they are not prepared for. This error may occur 4137when the matrix is specified to be real and the routines are not compiled 4138for real matrices, or when the matrix is specified to be complex and the 4139routines are not compiled to handle complex matrices. 4140 4141spFATAL - 2 4142 4143Not an error flag, but rather the dividing line between fatal errors and 4144warnings. 4145 4146 4147 4148 4149 4150 4151 4152 4153 June 23, 1988 4154 4155 4156 4157 4158 4159 - 62 - 4160 4161 41627.2: Data Structures 4163 4164 There is only one data structure that may need to be imported from 4165Sparse by the user. This data structure is used to hold pointers to four 4166related elements in matrix. It is used in conjunction with the routines 4167 spGetAdmittance() 4168 spGetOnes() 4169 spGetQuad() 4170 4171spGetAdmittance(), spGetOnes(), and spGetQuad() stuff the structure which 4172is later used by the spADD QUAD() macros. It is also possible for the user 4173to collect four pointers returned by spGetElement() and stuff them into the 4174template. The spADD QUAD() macros add a value into Element1 and Element2 4175and subtract the value from Element3 and Element4. The structure is: 4176 4177 4178struct spTemplate 4179{ spREAL *Element1; 4180 spREAL *Element2; 4181 spREAL *Element3Negated; 4182 spREAL *Element4Negated; 4183}; 4184 4185 4186 4187 4188 4189 4190 4191 4192 4193 4194 4195 4196 4197 4198 4199 4200 4201 4202 4203 4204 4205 4206 4207 4208 4209 4210 4211 4212 4213 4214 4215 4216 4217 4218 4219 4220 June 23, 1988 4221 4222 4223 4224 4225 4226 - 63 - 4227 4228 42298: FORTRAN COMPATIBILITY 4230 4231 The Sparse1.3 package contains routines that interface to a calling 4232program written in FORTRAN. Almost every externally available Sparse1.3 4233routine has a counterpart defined with the same name except that the `sp' 4234prefix is changed to `sf'. The spADD ELEMENT() and spADD QUAD() macros are 4235also replaced with the sfAdd1() and sfAdd4() functions. 4236 4237 Any interface between two languages is going to have portibility prob- 4238lems, this one is no exception. To ease porting the FORTRAN interface file 4239to different operating systems, the names of the interface functions can be 4240easily redefined (search for `Routine Renaming' in spFortran.c). When 4241interfacing to a FORTRAN program, the FORTRAN option should be set to YES 4242and the ARRAY OFFSET option should be set to NO (see spConfig.h). For 4243details on the return value and argument list of a particular interface 4244routine, see the file spFortran.c. 4245 4246 A simple example of a FORTRAN program that calls Sparse follows. 4247 4248 4249 4250 4251 4252 4253 4254 4255 4256 4257 4258 4259 4260 4261 4262 4263 4264 4265 4266 4267 4268 4269 4270 4271 4272 4273 4274 4275 4276 4277 4278 4279 4280 4281 4282 4283 4284 4285 4286 June 23, 1988 4287 4288 4289 4290 4291 4292 - 64 - 4293 4294 4295Example: 4296 integer matrix, error, sfCreate, sfGetElement, spFactor 4297 integer element(10) 4298 double precision rhs(4), solution(4) 4299 c 4300 c create matrix 4301 matrix = sfCreate(4,0,error) 4302 c 4303 c reserve elements 4304 element(1) = sfGetElement(matrix,1,1) 4305 element(2) = sfGetElement(matrix,1,2) 4306 element(3) = sfGetElement(matrix,2,1) 4307 element(4) = sfGetElement(matrix,2,2) 4308 element(5) = sfGetElement(matrix,2,3) 4309 element(6) = sfGetElement(matrix,3,2) 4310 element(7) = sfGetElement(matrix,3,3) 4311 element(8) = sfGetElement(matrix,3,4) 4312 element(9) = sfGetElement(matrix,4,3) 4313 element(10) = sfGetElement(matrix,4,4) 4314 c 4315 c clear matrix 4316 call sfClear(matrix) 4317 c 4318 c load matrix 4319 call sfAdd1Real(element(1), 2d0) 4320 call sfAdd1Real(element(2), -1d0) 4321 call sfAdd1Real(element(3), -1d0) 4322 call sfAdd1Real(element(4), 3d0) 4323 call sfAdd1Real(element(5), -1d0) 4324 call sfAdd1Real(element(6), -1d0) 4325 call sfAdd1Real(element(7), 3d0) 4326 call sfAdd1Real(element(8), -1d0) 4327 call sfAdd1Real(element(9), -1d0) 4328 call sfAdd1Real(element(10), 3d0) 4329 call sfprint(matrix, .false., .false., .true.) 4330 rhs(1) = 34d0 4331 rhs(2) = 0d0 4332 rhs(3) = 0d0 4333 rhs(4) = 0d0 4334 c 4335 c factor matrix 4336 error = sfFactor(matrix) 4337 c 4338 c solve matrix 4339 call sfSolve(matrix, rhs, solution) 4340 write (6, 10) solution(1), solution(2), solution(3), solution(4) 4341 10 format (f 10.2) 4342 end 4343 4344 4345 4346 4347 4348 4349 4350 4351 4352 June 23, 1988 4353 4354 4355 4356 4357 4358 - 65 - 4359 4360 43619: SPARSE TEST PROGRAM 4362 4363 The Sparse package includes a test program that is able to read matrix 4364equations from text files and print their solution along with matrix 4365statistics and timing information. The program can also generate files 4366containing stripped versions of the unfactored and factored matrix suitable 4367for plotting using standard plotting programs, such as the UNIX graph and 4368plot commands. 4369 4370The Sparse test program is invoked using the following syntax. 4371 4372 sparse [options] [file1] [file2] ... 4373 4374 Options: 4375 -s Print solution only. 4376 -r x Use x as relative threshold. 4377 -a x Use x as absolute threshold. 4378 -n n Print first n terms of solution vector. 4379 -i n Repeat build/factor/solve n times for better 4380 timing results. 4381 -b n Use column n of matrix as right-hand side 4382 vector. 4383 -p Create plot files ``filename.bef'' and 4384 ``filename.aft''. 4385 -c Use complete (as opposed to diagonal) pivot- 4386 ing. 4387 -x Treat real matrix as complex with imaginary 4388 part zero. 4389 -t Solve transposed system. 4390 -u Print usage message. 4391 4392 4393The presence of certain options is dependent on whether the appropriate 4394Sparse option has been enabled. 4395 4396If no input files are specified, sparse reads from the standard input. The 4397syntax of the input file is as follows. The matrix begins with one line of 4398arbitrary text that acts as the label, followed by a line with the integer 4399size of the matrix and either the real or complex keywords. After the 4400header is an arbitrary number of lines that describe the structural 4401nonzeros in the matrix. These lines have the form row column data, where 4402row and column are integers and data is either one real number for real 4403matrices or a real/imaginary pair of numbers for complex matrices. Only 4404one structural nonzero is described per line and the section ends when 4405either row or column are zero. Following the matrix, an optional right- 4406hand side vector can be described. The vector is given one element per 4407line, the number of element must equal the size of the matrix. Only one 4408matrix and one vector are allowed per file, and the vector, if given, must 4409follow the matrix. 4410 4411 4412 4413 4414 4415 4416 4417 4418 June 23, 1988 4419 4420 4421 4422 4423 4424 - 66 - 4425 4426 4427Example: 4428 mat0 - Simple matrix. 4429 4 real 4430 1 1 2.0 4431 1 2 -1.0 4432 2 1 -1.0 4433 2 2 3.0 4434 2 3 -1.0 4435 3 2 -1.0 4436 3 3 3.0 4437 3 4 -1.0 4438 4 3 -1.0 4439 4 4 3.0 4440 0 0 0.0 4441 34.0 4442 0.0 4443 0.0 4444 0.0 4445 4446 4447 4448 4449 4450 4451 4452 4453 4454 4455 4456 4457 4458 4459 4460 4461 4462 4463 4464 4465 4466 4467 4468 4469 4470 4471 4472 4473 4474 4475 4476 4477 4478 4479 4480 4481 4482 4483 4484 June 23, 1988 4485 4486 4487 4488 4489 4490 - 67 - 4491 4492 449310: SPARSE FILES 4494 4495 The following is a list of the files contained in the Sparse package 4496and a brief description of their contents. Of the files, only spConfig.h 4497is expected to be modified by the user and only spMatrix.h need be imported 4498into the program that calls Sparse. 4499 4500 4501spAlloc.c 4502 4503This file contains the routines for allocating and deallocating objects as- 4504sociated with the matrices, including the matrices themselves. 4505 4506o User accessible functions contained in this module: 4507 spCreate() 4508 spDestroy() 4509 spError() 4510 spWhereSingular() 4511 spGetSize() 4512 spSetReal() 4513 spSetComplex() 4514 spFillinCount() 4515 spElementCount() 4516 4517 4518spBuild.c 4519 4520This file contains the routines for clearing and loading the matrix. 4521 4522o User accessible functions contained in this module: 4523 spClear() 4524 spGetAdmittance() 4525 spGetElement() 4526 spGetInitInfo() 4527 spGetOnes() 4528 spGetQuad() 4529 spInitialize() 4530 spInstallInitInfo() 4531 4532 4533 4534 4535 4536 4537 4538 4539 4540 4541 4542 4543 4544 4545 4546 4547 4548 4549 4550 June 23, 1988 4551 4552 4553 4554 4555 4556 - 68 - 4557 4558 4559 4560spCompat.c 4561 4562This file contains the routines for making Sparse1.3 upward compatible from 4563Sparse1.2. These routines are not suggested for use in new software. 4564These routines will not be available in future versions of Sparse. 4565 4566o User accessible functions contained in this module: 4567 AddAdmittanceToMatrix() 4568 AddComplexElementToMatrix() 4569 AddComplexQuadElementToMatrix() 4570 AddElementToMatrix() 4571 AddImagElementToMatrix() 4572 AddImagQuadElementToMatrix() 4573 AddOnesToMatrix() 4574 AddQuadToMatrix() 4575 AddRealElementToMatrix() 4576 AddRealQuadElementToMatrix() 4577 CleanMatrix() 4578 ClearMatrix() 4579 ClearMatrixError() 4580 CreateMatrix() 4581 DecomposeMatrix() 4582 DeleteRowAndColFromMatrix() 4583 DestroyMatrix() 4584 Determinant() 4585 GetMatrixSize() 4586 MatrixElementCount() 4587 MatrixError() 4588 MatrixFillinCount() 4589 MatrixRoundoffError() 4590 MultiplyMatrix() 4591 OrderAndDecomposeMatrix() 4592 OutputMatrixToFile() 4593 PreorderForModifiedNodal() 4594 PrintMatrix() 4595 ScaleMatrix() 4596 SetMatrixComplex() 4597 SetMatrixReal() 4598 SolveMatrix() 4599 SolveTransposedMatrix() 4600 4601 4602spConfig.h 4603 4604This file contains the options that are used to customize the package. For 4605example, it is possible to specify whether only real or complex systems of 4606equations are to be solved. Also included in this file are the various 4607constants used by the Sparse package, such as the amount of memory initial- 4608ly allocated for each matrix and the largest real number represented by the 4609machine. The user is expected to modify this file to maximize the perfor- 4610mance of the routines with his/her matrices. 4611 4612 4613 4614 4615 4616 June 23, 1988 4617 4618 4619 4620 4621 4622 - 69 - 4623 4624 4625 4626spDefs.h 4627 4628This module contains common data structure definitions and macros for the 4629sparse matrix routines. These definitions are meant to remain hidden from 4630the program that calls the sparse matrix routines. 4631 4632 4633spDoc 4634 4635This reference manual. spDoc contains the manual in a form that is read- 4636able on-line and spDoc.ms contains the manual in a form that is suitable 4637for input into the text formatting program troff using the -ms macros. 4638 4639 4640spFactor.c 4641 4642This file contains the routines for factoring matrices into LU form. 4643 4644o User accessible functions contained in this module: 4645 spFactor() 4646 spOrderAndFactor() 4647 spPartition() 4648 4649 4650 4651 4652 4653 4654 4655 4656 4657 4658 4659 4660 4661 4662 4663 4664 4665 4666 4667 4668 4669 4670 4671 4672 4673 4674 4675 4676 4677 4678 4679 4680 4681 4682 June 23, 1988 4683 4684 4685 4686 4687 4688 - 70 - 4689 4690 4691 4692spFortran.c 4693 4694This file contains the routines for interfacing Sparse1.3 to a program 4695written in FORTRAN. The function and argument lists of the routines in 4696this file are almost identical to their C equivalents except that they are 4697suitable for calling from a FORTRAN program. The names of these routines 4698use the `sf' prefix to distinguish them from their C counterparts. 4699 4700o User accessible functions contained in this module: 4701 sfAdd1Complex() 4702 sfAdd1Imag() 4703 sfAdd1Real() 4704 sfAdd4Complex() 4705 sfAdd4Imag() 4706 sfAdd4Real() 4707 sfClear() 4708 sfCondition() 4709 sfCreate() 4710 sfDeleteRowAndCol() 4711 sfDestroy() 4712 sfDeterminant() 4713 sfElementCount() 4714 sfError() 4715 sfFactor() 4716 sfFileMatrix() 4717 sfFileStats() 4718 sfFileVector() 4719 sfFillinCount() 4720 sfGetAdmittance() 4721 sfGetElement() 4722 sfGetOnes() 4723 sfGetQuad() 4724 sfGetSize() 4725 sfLargestElement() 4726 sfMNA Preorder() 4727 sfMultTransposed() 4728 sfMultiply() 4729 sfNorm() 4730 sfOrderAndFactor() 4731 sfPartition() 4732 sfPrint() 4733 sfPseudoCondition() 4734 sfRoundoff() 4735 sfScale() 4736 sfSetComplex() 4737 sfSetReal() 4738 sfSolve() 4739 sfSolveTransposed() 4740 sfStripFills() 4741 sfWhereSingular() 4742 4743 4744 4745 4746 4747 4748 June 23, 1988 4749 4750 4751 4752 4753 4754 - 71 - 4755 4756 4757 4758spMatrix.h 4759 4760This file contains definitions that are useful to the calling program. In 4761particular, this file contains error keyword definitions, some macro func- 4762tions that are used to quickly enter data into the matrix, the definition 4763of a data structure that acts as a template for entering admittances into 4764the matrix, and the type declarations of the various Sparse functions. 4765 4766 4767spOutput.c 4768 4769This file contains the output-to-file and output-to-screen routines for the 4770matrix package. They are capable of outputting the matrix in either a form 4771readable by people or a form readable by the Sparse test program. 4772 4773o User accessible functions contained in this module: 4774 spFileMatrix() 4775 spFileStats() 4776 spFileVector() 4777 spPrint() 4778 4779 4780spRevision 4781 4782The history of updates for the program. This file also includes ordering 4783information for the Sparse package. 4784 4785 4786spSolve.c 4787 4788This module contains the forward and backward substitution routines. 4789 4790o User accessible functions contained in this module: 4791 spSolve() 4792 spSolveTransposed() 4793 4794 4795spTest.c 4796 4797This module contains a test program for the sparse matrix routines. It is 4798able to read matrices from files and solve them. Because of the large 4799number of options and capabilities built into Sparse, it is impossible to 4800have one test routine thoroughly exercise Sparse. Thus, emphasis is on ex- 4801ercising as many capabilities as is reasonable while also providing a use- 4802ful tool. 4803 4804 4805 4806 4807 4808 4809 4810 4811 4812 4813 4814 June 23, 1988 4815 4816 4817 4818 4819 4820 - 72 - 4821 4822 4823 4824spUtil.c 4825 4826This module contains various optional utility routines. 4827 4828o User accessible functions contained in this module: 4829 spCondition() 4830 spDeleteRowAndCol() 4831 spDeterminant() 4832 spLargestElement() 4833 spMNA Preorder() 4834 spMultiply() 4835 spMultTransposed() 4836 spNorm() 4837 spPseudoCondition() 4838 spRoundoff() 4839 spScale() 4840 spStripFills() 4841 4842 4843Makefile 4844 4845This file is used in conjunction with the UNIX program make to compile the 4846matrix routines and their test program. 4847 4848 4849make.com 4850 4851This file is used to automatically compile Sparse under the VMS operating 4852system. It needs to modified slightly before being used, see the installa- 4853tion notes. 4854 4855 4856 4857 4858 4859 4860 4861 4862 4863 4864 4865 4866 4867 4868 4869 4870 4871 4872 4873 4874 4875 4876 4877 4878 4879 4880 June 23, 1988 4881 4882 4883 4884 4885 4886 - 73 - 4887 4888 4889REFERENCES 4890 4891[duff86] I. S. Duff, A. M. Erisman, J. K. Reid. Direct Methods for 4892 Sparse Matrices. Oxford University Press, 1986. 4893 4894[golub86] G. H. Golub, C. F. V. Van Loan. Matrix Computations. The 4895 Johns Hopkins University Press, 1983. 4896 4897[kundert86] Kenneth S. Kundert. Sparse matrix techniques. In Circuit 4898 Analysis, Simulation and Design, Albert Ruehli (editor). 4899 North-Holland, 1986. 4900 4901[strang80] Gilbert Strang. Linear Algebra and Its Applications. 4902 Academic Press, 1980. 4903 4904 4905Acknowledgements 4906 4907 We would like to acknowledge and thank the those people that contri- 4908buted ideas that were incorporated into Sparse. In particular, Jacob 4909White, Kartikeya Mayaram, Don Webber, Tom Quarles, Howard Ko and Beresford 4910Parlett. 4911 4912 4913 4914 4915 4916 4917 4918 4919 4920 4921 4922 4923 4924 4925 4926 4927 4928 4929 4930 4931 4932 4933 4934 4935 4936 4937 4938 4939 4940 4941 4942 4943 4944 4945 4946 June 23, 1988 4947 4948 4949 4950 4951 4952 4953 4954 4955 4956 4957 4958 Table of Contents 4959 4960 4961 4962 49631: Introduction ..................................................... 1 4964 4965 1.1: Features of Sparse1.3 .................................. 1 4966 4967 1.2: Enhancements of Sparse1.3 over Sparse1.2 ............... 2 4968 4969 1.3: Copyright Information .................................. 3 4970 49712: Primer ........................................................... 4 4972 4973 2.1: Solving Matrix Equations ............................... 4 4974 4975 2.2: Error Control .......................................... 5 4976 4977 2.3: Building the Matrix .................................... 6 4978 4979 2.4: Initializing the Matrix ................................ 7 4980 4981 2.5: Indices ................................................ 8 4982 4983 2.6: Configuring Sparse ..................................... 9 4984 49853: Introduction to the Sparse Routines .............................. 10 4986 4987 3.1: Creating the Matrix .................................... 10 4988 4989 3.2: Building the Matrix .................................... 10 4990 4991 3.3: Clearing the Matrix .................................... 10 4992 4993 3.4: Placing Data in the Matrix ............................. 11 4994 4995 3.5: Influencing the Factorization .......................... 11 4996 4997 3.6: Factoring the Matrix ................................... 11 4998 4999 3.7: Solving the Matrix Equation ............................ 12 5000 5001 3.8: Numerical Error Estimation ............................. 12 5002 5003 3.9: Matrix Operations ...................................... 13 5004 5005 3.10: Matrix Statistics and Documentation ................... 13 5006 50074: Routines ......................................................... 15 5008 5009 5010 5011 5012 June 23, 1988 5013 5014 5015 5016 5017 5018 5019 5020 5021 4.1: spClear() .............................................. 15 5022 5023 4.2: spCondition() .......................................... 16 5024 5025 4.3: spCreate() ............................................. 17 5026 5027 4.4: spDeleteRowAndCol() .................................... 18 5028 5029 4.5: spDestroy() ............................................ 18 5030 5031 4.6: spDeterminant() ........................................ 19 5032 5033 4.7: spElementCount() ....................................... 20 5034 5035 4.8: spError() .............................................. 20 5036 5037 4.9: spFactor() ............................................. 21 5038 5039 4.10: spFileMatrix() ........................................ 22 5040 5041 4.11: spFileStats() ......................................... 23 5042 5043 4.12: spFileVector() ........................................ 24 5044 5045 4.13: spFillinCount() ....................................... 25 5046 5047 4.14: spGetAdmittance() ..................................... 26 5048 5049 4.15: spGetElement() ........................................ 27 5050 5051 4.16: spGetInitInfo() ....................................... 28 5052 5053 4.17: spGetOnes() ........................................... 30 5054 5055 4.18: spGetQuad() ........................................... 32 5056 5057 4.19: spGetSize() ........................................... 32 5058 5059 4.20: spInitialize() ........................................ 34 5060 5061 4.21: spInstallInitInfo() ................................... 34 5062 5063 4.22: spLargestElement() .................................... 35 5064 5065 4.23: spMNA Preorder() ...................................... 36 5066 5067 4.24: spMultiply() .......................................... 37 5068 5069 4.25: spMultTransposed() .................................... 38 5070 5071 4.26: spNorm() .............................................. 39 5072 5073 4.27: spOrderAndFactor() .................................... 39 5074 5075 5076 5077 5078 June 23, 1988 5079 5080 5081 5082 5083 5084 5085 5086 5087 4.28: spPartition() ......................................... 42 5088 5089 4.29: spPrint() ............................................. 42 5090 5091 4.30: spPseudoCondition() ................................... 43 5092 5093 4.31: spRoundoff() .......................................... 44 5094 5095 4.32: spScale() ............................................. 45 5096 5097 4.33: spSetComplex() ........................................ 46 5098 5099 4.34: spSetReal() ........................................... 46 5100 5101 4.35: spSolve() ............................................. 47 5102 5103 4.36: spSolveTransposed() ................................... 48 5104 5105 4.37: spStripFills() ........................................ 49 5106 5107 4.38: spWhereSingular() ..................................... 49 5108 51095: Macro Functions .................................................. 50 5110 5111 5.1: spADD REAL ELEMENT() ................................... 50 5112 5113 5.2: spADD IMAG ELEMENT() ................................... 50 5114 5115 5.3: spADD COMPLEX ELEMENT() ................................ 51 5116 5117 5.4: spADD REAL QUAD() ...................................... 51 5118 5119 5.5: spADD IMAG QUAD() ...................................... 52 5120 5121 5.6: spADD COMPLEX QUAD() ................................... 52 5122 51236: Configuring Sparse ............................................... 53 5124 5125 6.1: Sparse Options ......................................... 53 5126 5127 6.2: Sparse Constants ....................................... 58 5128 5129 6.3: Machine Constants ...................................... 60 5130 51317: Exports .......................................................... 61 5132 5133 7.1: Error Codes ............................................ 61 5134 5135 7.2: Data Structures ........................................ 62 5136 51378: FORTRAN Compatibility ............................................ 63 5138 51399: Sparse Test Program .............................................. 65 5140 5141 5142 5143 5144 June 23, 1988 5145 5146 5147 5148 5149 5150 5151 5152 515310: Sparse Files .................................................... 67 5154 5155References ........................................................... 73 5156 5157 5158 5159 5160 5161 5162 5163 5164 5165 5166 5167 5168 5169 5170 5171 5172 5173 5174 5175 5176 5177 5178 5179 5180 5181 5182 5183 5184 5185 5186 5187 5188 5189 5190 5191 5192 5193 5194 5195 5196 5197 5198 5199 5200 5201 5202 5203 5204 5205 5206 5207 5208 5209 5210 June 23, 1988 5211 5212 5213