home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol244 / lpak11-c.lqr / SGEFA.C < prev    next >
Encoding:
Text File  |  1986-02-11  |  5.6 KB  |  123 lines
  1. /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */
  2.  
  3. SGEFA ( A, lda, n, IPvt, InfoCode )
  4.    float A[][ldx] ;
  5.    int lda, n ;
  6.    int IPvt[] ;
  7.    int *InfoCode ;
  8.  
  9.      /* Program:                                            */
  10.      /*                                                     */
  11.      /*    SGEFA - Factors Real Single Precision Matrix     */
  12.      /*            Using Gaussian Elimination.              */
  13.      /*                                                     */
  14.      /* Version:                 DATE:                      */
  15.      /*                                                     */
  16.      /*     C                       02/22/85                */
  17.      /*     FORTRAN                 08/14/78                */
  18.      /*                                                     */
  19.      /* Description:                                        */
  20.      /*                                                     */
  21.      /*     SGEFA is usually called by SGECO, but can be    */
  22.      /*     called directly with a saving in time if RCOND  */
  23.      /*     is not needed.                                  */
  24.      /*                                                     */
  25.      /*     (time for SGECO) = (1 + 9/n)*(time for SGEFA)   */
  26.      /*                                                     */
  27.      /*     On entry;                                       */
  28.      /*                                                     */
  29.      /*        A     - Matrix to be factored.               */
  30.      /*        lda   - Leading dimension of matrix A.       */
  31.      /*        n     - Order of matrix A.                   */
  32.      /*                                                     */
  33.      /*     On exit;                                        */
  34.      /*                                                     */
  35.      /*        A     - Upper triangular matrix and multi-   */
  36.      /*                pliers used to obtain it. Factoriza- */
  37.      /*                tion can be written;                 */
  38.      /*                   A = L * U                         */
  39.      /*                where L is a product of permutation  */
  40.      /*                and unit lower triangular matrices   */
  41.      /*                and U is upper triangular.           */
  42.      /*                                                     */
  43.      /*        IPvt  - Pivot index vector.                  */
  44.      /*                                                     */
  45.      /*        InfoCode -                                   */
  46.      /*            = 0 Normal value.                        */
  47.      /*            = k If  u[k,k] = 0.0. This is not an     */
  48.      /*                error condition for this procedure,  */
  49.      /*                but indicates that SGESL or SGEDI    */
  50.      /*                will divide by zero if called.  Use  */
  51.      /*                RCOND from SGECO for reliable indi-  */
  52.      /*                cation of singularity.               */
  53.      /*                                                     */
  54.      /* Subprograms:                                        */
  55.      /*                                                     */
  56.      /*    SAXPY       from BLAS                            */
  57.      /*    SSCAL       from BLAS                            */
  58.      /*    ISAMAX      from BLAS                            */
  59.      /*                                                     */
  60.      /* Authors:                                            */
  61.      /*                                                     */
  62.      /*    Cleve Moler                                      */
  63.      /*    University of New Mexico                         */
  64.      /*    Argonne National Laboratories                    */
  65.      /*                                                     */
  66.      /*    C and Pascal translations;                       */
  67.      /*                                                     */
  68.      /*    Adam Fritz                                       */
  69.      /*    133 Main Street                                  */
  70.      /*    Afton, New York 13730                            */
  71.  
  72. {
  73.    int ISAMAX() ;
  74.    int j, k, l ;
  75.    int kp1, nm1 ;
  76.    float t ;
  77.                                 /* Gaussian Elimination with Partial Pivoting */
  78.    *InfoCode = 0 ;
  79.    nm1 = n - 1 ;
  80.    if ( nm1 > 0 )
  81.    {
  82.       for ( k = 0; k < nm1; k++ )
  83.       {
  84.          kp1 = k + 1 ;
  85.                                 /* Find l = Pivot Index */
  86.          l = ISAMAX(n-k, &A[k][k], lda) + k ;
  87.          IPvt[k] = l ;
  88.                                 /* Zero Pivot Implies Column Triangularized */
  89.          if ( A[l][k] != 0.0 )
  90.          {
  91.                                 /* Interchange If Necessary */
  92.             if ( l != k )
  93.             {
  94.                t = A[l][k] ;
  95.                A[l][k] = A[k][k] ;
  96.                A[k][k] = t ;
  97.             }
  98.                                 /* Compute Multipliers */
  99.             t = -1.0/A[k][k] ;
  100.             SSCAL(n-k-1, t, &A[k+1][k], lda) ;
  101.                                 /* Row Elimination with Column Indexing */
  102.             for ( j = kp1; j < n; j++ )
  103.             {
  104.                t = A[l][j] ;
  105.                if ( l != k )
  106.                {
  107.                   A[l][j] = A[k][j] ;
  108.                   A[k][j] = t ;
  109.                }
  110.                SAXPY(n-k-1, t, &A[k+1][k], lda, &A[k+1][j], lda) ;
  111.             }
  112.          }
  113.          else
  114.             *InfoCode = k ;
  115.       }
  116.    }
  117.    IPvt[nm1] = nm1 ;
  118.    if ( A[nm1][nm1] == 0.0 )
  119.       *InfoCode = nm1 ;
  120. }
  121.  
  122. /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */
  123.