Return to linpack.c CVS log

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

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

1.2     ! albertel    1: /* functions to support the random number genrator
        !             2:    Copyright (C) 1992-2000 Michigan State University
        !             3: 
        !             4:    The CAPA system is free software; you can redistribute it and/or
        !             5:    modify it under the terms of the GNU Library General Public License as
        !             6:    published by the Free Software Foundation; either version 2 of the
        !             7:    License, or (at your option) any later version.
        !             8: 
        !             9:    The CAPA system is distributed in the hope that it will be useful,
        !            10:    but WITHOUT ANY WARRANTY; without even the implied warranty of
        !            11:    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
        !            12:    Library General Public License for more details.
        !            13: 
        !            14:    You should have received a copy of the GNU Library General Public
        !            15:    License along with the CAPA system; see the file COPYING.  If not,
        !            16:    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
        !            17:    Boston, MA 02111-1307, USA.  */
        !            18: 
1.1       albertel   19: #include <math.h>
                     20: float sdot(long n,float *sx,long incx,float *sy,long incy)
                     21: {
                     22: static long i,ix,iy,m,mp1;
                     23: static float sdot,stemp;
                     24:     stemp = sdot = 0.0;
                     25:     if(n <= 0) return sdot;
                     26:     if(incx == 1 && incy == 1) goto S20;
                     27:     ix = iy = 1;
                     28:     if(incx < 0) ix = (-n+1)*incx+1;
                     29:     if(incy < 0) iy = (-n+1)*incy+1;
                     30:     for(i=1; i<=n; i++) {
                     31:         stemp += (*(sx+ix-1)**(sy+iy-1));
                     32:         ix += incx;
                     33:         iy += incy;
                     34:     }
                     35:     sdot = stemp;
                     36:     return sdot;
                     37: S20:
                     38:     m = n % 5L;
                     39:     if(m == 0) goto S40;
                     40:     for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
                     41:     if(n < 5) goto S60;
                     42: S40:
                     43:     mp1 = m+1;
                     44:     for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
                     45:       +1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
                     46: S60:
                     47:     sdot = stemp;
                     48:     return sdot;
                     49: }
                     50: void spofa(float *a,long lda,long n,long *info)
                     51: /*
                     52:      SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
                     53:      SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
                     54:      DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
                     55:      (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
                     56:      ON ENTRY
                     57:         A       REAL(LDA, N)
                     58:                 THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
                     59:                 DIAGONAL AND UPPER TRIANGLE ARE USED.
                     60:         LDA     INTEGER
                     61:                 THE LEADING DIMENSION OF THE ARRAY  A .
                     62:         N       INTEGER
                     63:                 THE ORDER OF THE MATRIX  A .
                     64:      ON RETURN
                     65:         A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
                     66:                 WHERE  TRANS(R)  IS THE TRANSPOSE.
                     67:                 THE STRICT LOWER TRIANGLE IS UNALTERED.
                     68:                 IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
                     69:         INFO    INTEGER
                     70:                 = 0  FOR NORMAL RETURN.
                     71:                 = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
                     72:                      OF ORDER  K  IS NOT POSITIVE DEFINITE.
                     73:      LINPACK.  THIS VERSION DATED 08/14/78 .
                     74:      CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
                     75:      SUBROUTINES AND FUNCTIONS
                     76:      BLAS SDOT
                     77:      FORTRAN SQRT
                     78:      INTERNAL VARIABLES
                     79: */
                     80: {
                     81: extern float sdot(long n,float *sx,long incx,float *sy,long incy);
                     82: static long j,jm1,k;
                     83: static float t,s;
                     84: /*
                     85:      BEGIN BLOCK WITH ...EXITS TO 40
                     86: */
                     87:     for(j=1; j<=n; j++) {
                     88:         *info = j;
                     89:         s = 0.0;
                     90:         jm1 = j-1;
                     91:         if(jm1 < 1) goto S20;
                     92:         for(k=0; k<jm1; k++) {
                     93:             t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
                     94:             t /=  *(a+k+k*lda);
                     95:             *(a+k+(j-1)*lda) = t;
                     96:             s += (t*t);
                     97:         }
                     98: S20:
                     99:         s = *(a+j-1+(j-1)*lda)-s;
                    100: /*
                    101:      ......EXIT
                    102: */
                    103:         if(s <= 0.0) goto S40;
                    104:         *(a+j-1+(j-1)*lda) = sqrt(s);
                    105:     }
                    106:     *info = 0;
                    107: S40:
                    108:     return;
                    109: }