/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:54 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sgbfa.h" #include void /*FUNCTION*/ sgbfa( float *abd, long lda, long n, long ml, long mu, long ipvt[], long *info) { #define ABD(I_,J_) (*(abd+(I_)*(lda)+(J_))) long int i, i0, j, j0, j1, ju, jz, k, kp1, l, lm, m, mm, nm1; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ipvt = &ipvt[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-11-04 SGBFA Krogh Fixes for F77 and conversion to single *--S replaces "?": ?GBFA, ?GBCO, ?AXPY, ?SCAL, I?AMAX * IMPLICIT NONE ****BEGIN PROLOGUE SGBFA ****PURPOSE Factor a band matrix using Gaussian elimination. ****LIBRARY SLATEC (LINPACK) ****CATEGORY D2A2 ****TYPE DOUBLE PRECISION (SGBFA-S, SGBFA-D, CGBFA-C) ****KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION ****AUTHOR Moler, C. B., (U. of New Mexico) ****DESCRIPTION * * SGBFA factors a double precision band matrix by elimination. * * SGBFA is usually called by SGBCO, but it can be called * directly with a saving in time if RCOND is not needed. * * On Entry * * ABD DOUBLE PRECISION(LDA, N) * contains the matrix in band storage. The columns * of the matrix are stored in the columns of ABD and * the diagonals of the matrix are stored in rows * ML+1 through 2*ML+MU+1 of ABD . * See the comments below for details. * * LDA INTEGER * the leading dimension of the array ABD . * LDA must be .GE. 2*ML + MU + 1 . * * N INTEGER * the order of the original matrix. * * ML INTEGER * number of diagonals below the main diagonal. * 0 .LE. ML .LT. N . * * MU INTEGER * number of diagonals above the main diagonal. * 0 .LE. MU .LT. N . * More efficient if ML .LE. MU . * On Return * * ABD an upper triangular matrix in band storage and * the multipliers which were used to obtain it. * The factorization can be written A = L*U where * L is a product of permutation and unit lower * triangular matrices and U is upper triangular. * * IPVT INTEGER(N) * an integer vector of pivot indices. * * INFO INTEGER * = 0 normal value. * = K if U(K,K) .EQ. 0.0 . This is not an error * condition for this subroutine, but it does * indicate that DGBSL will divide by zero if * called. Use RCOND in SGBCO for a reliable * indication of singularity. * * Band Storage * * If A is a band matrix, the following program segment * will set up the input. * * ML = (band width below the diagonal) * MU = (band width above the diagonal) * M = ML + MU + 1 * DO 20 J = 1, N * I1 = MAX(1, J-MU) * I2 = MIN(N, J+ML) * DO 10 I = I1, I2 * K = I - J + M * ABD(K,J) = A(I,J) * 10 CONTINUE * 20 CONTINUE * * This uses rows ML+1 through 2*ML+MU+1 of ABD . * In addition, the first ML rows in ABD are used for * elements generated during the triangularization. * The total number of rows needed in ABD is 2*ML+MU+1 . * The ML+MU by ML+MU upper left triangle and the * ML by ML lower right triangle are not referenced. * ****REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. * Stewart, LINPACK Users' Guide, SIAM, 1979. ****ROUTINES CALLED SAXPY, SSCAL, ISAMAX ****REVISION HISTORY (YYMMDD) * 780814 DATE WRITTEN * 890531 Changed all specific intrinsics to generic. (WRB) * 890831 Modified array declarations. (WRB) * 890831 REVISION DATE from Version 3.2 * 891214 Prologue converted to Version 4.0 format. (BAB) * 900326 Removed duplicate information from DESCRIPTION section. * (WRB) * 920501 Reformatted the REFERENCES section. (WRB) ****END PROLOGUE SGBFA */ /****FIRST EXECUTABLE STATEMENT SGBFA */ m = ml + mu + 1; *info = 0; /* ZERO INITIAL FILL-IN COLUMNS * */ j0 = mu + 2; j1 = min( n, m ) - 1; if (j1 < j0) goto L_30; for (jz = j0; jz <= j1; jz++) { i0 = m + 1 - jz; for (i = i0; i <= ml; i++) { ABD(jz - 1,i - 1) = 0.0e0; } } L_30: ; jz = j1; ju = 0; /* GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING * */ nm1 = n - 1; if (nm1 < 1) goto L_130; for (k = 1; k <= nm1; k++) { kp1 = k + 1; /* ZERO NEXT FILL-IN COLUMN * */ jz += 1; if (jz > n) goto L_50; if (ml < 1) goto L_50; for (i = 1; i <= ml; i++) { ABD(jz - 1,i - 1) = 0.0e0; } L_50: ; /* FIND L = PIVOT INDEX * */ lm = min( ml, n - k ); l = isamax( lm + 1, &ABD(k - 1,m - 1), 1 ) + m - 1; Ipvt[k] = l + k - m; /* ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED * */ if (ABD(k - 1,l - 1) == 0.0e0) goto L_100; /* INTERCHANGE IF NECESSARY * */ if (l == m) goto L_60; t = ABD(k - 1,l - 1); ABD(k - 1,l - 1) = ABD(k - 1,m - 1); ABD(k - 1,m - 1) = t; L_60: ; /* COMPUTE MULTIPLIERS * */ t = -1.0e0/ABD(k - 1,m - 1); sscal( lm, t, &ABD(k - 1,m), 1 ); /* ROW ELIMINATION WITH COLUMN INDEXING * */ ju = min( max( ju, mu + Ipvt[k] ), n ); mm = m; if (ju < kp1) goto L_90; for (j = kp1; j <= ju; j++) { l -= 1; mm -= 1; t = ABD(j - 1,l - 1); if (l == mm) goto L_70; ABD(j - 1,l - 1) = ABD(j - 1,mm - 1); ABD(j - 1,mm - 1) = t; L_70: ; saxpy( lm, t, &ABD(k - 1,m), 1, &ABD(j - 1,mm), 1 ); } L_90: ; goto L_110; L_100: ; *info = k; L_110: ; } L_130: ; Ipvt[n] = n; if (ABD(n - 1,m - 1) == 0.0e0) *info = n; return; #undef ABD } /* end of function */