Return to linpack.c CVS log

Up to [LON-CAPA] / capa / capa51 / pProj

Initial revision

    1: #include <math.h>
    2: float sdot(long n,float *sx,long incx,float *sy,long incy)
    3: {
    4: static long i,ix,iy,m,mp1;
    5: static float sdot,stemp;
    6:     stemp = sdot = 0.0;
    7:     if(n <= 0) return sdot;
    8:     if(incx == 1 && incy == 1) goto S20;
    9:     ix = iy = 1;
   10:     if(incx < 0) ix = (-n+1)*incx+1;
   11:     if(incy < 0) iy = (-n+1)*incy+1;
   12:     for(i=1; i<=n; i++) {
   13:         stemp += (*(sx+ix-1)**(sy+iy-1));
   14:         ix += incx;
   15:         iy += incy;
   16:     }
   17:     sdot = stemp;
   18:     return sdot;
   19: S20:
   20:     m = n % 5L;
   21:     if(m == 0) goto S40;
   22:     for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
   23:     if(n < 5) goto S60;
   24: S40:
   25:     mp1 = m+1;
   26:     for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
   27:       +1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
   28: S60:
   29:     sdot = stemp;
   30:     return sdot;
   31: }
   32: void spofa(float *a,long lda,long n,long *info)
   33: /*
   34:      SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
   35:      SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
   36:      DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
   37:      (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
   38:      ON ENTRY
   39:         A       REAL(LDA, N)
   40:                 THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
   41:                 DIAGONAL AND UPPER TRIANGLE ARE USED.
   42:         LDA     INTEGER
   43:                 THE LEADING DIMENSION OF THE ARRAY  A .
   44:         N       INTEGER
   45:                 THE ORDER OF THE MATRIX  A .
   46:      ON RETURN
   47:         A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
   48:                 WHERE  TRANS(R)  IS THE TRANSPOSE.
   49:                 THE STRICT LOWER TRIANGLE IS UNALTERED.
   50:                 IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
   51:         INFO    INTEGER
   52:                 = 0  FOR NORMAL RETURN.
   53:                 = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
   54:                      OF ORDER  K  IS NOT POSITIVE DEFINITE.
   55:      LINPACK.  THIS VERSION DATED 08/14/78 .
   56:      CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
   57:      SUBROUTINES AND FUNCTIONS
   58:      BLAS SDOT
   59:      FORTRAN SQRT
   60:      INTERNAL VARIABLES
   61: */
   62: {
   63: extern float sdot(long n,float *sx,long incx,float *sy,long incy);
   64: static long j,jm1,k;
   65: static float t,s;
   66: /*
   67:      BEGIN BLOCK WITH ...EXITS TO 40
   68: */
   69:     for(j=1; j<=n; j++) {
   70:         *info = j;
   71:         s = 0.0;
   72:         jm1 = j-1;
   73:         if(jm1 < 1) goto S20;
   74:         for(k=0; k<jm1; k++) {
   75:             t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
   76:             t /=  *(a+k+k*lda);
   77:             *(a+k+(j-1)*lda) = t;
   78:             s += (t*t);
   79:         }
   80: S20:
   81:         s = *(a+j-1+(j-1)*lda)-s;
   82: /*
   83:      ......EXIT
   84: */
   85:         if(s <= 0.0) goto S40;
   86:         *(a+j-1+(j-1)*lda) = sqrt(s);
   87:     }
   88:     *info = 0;
   89: S40:
   90:     return;
   91: }