Return to linpack.c CVS log

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

Annotation of capa/capa51/pProj/linpack.c, revision 1.1

1.1     ! albertel    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: }