Return to linpack.c CVS log

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

- updating GPL notices

    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: 
   19:    As a special exception, you have permission to link this program
   20:    with the TtH/TtM library and distribute executables, as long as you
   21:    follow the requirements of the GNU GPL in regard to all of the
   22:    software in the executable aside from TtH/TtM.
   23: */
   24: 
   25: #include <math.h>
   26: float sdot(long n,float *sx,long incx,float *sy,long incy)
   27: {
   28: static long i,ix,iy,m,mp1;
   29: static float sdot,stemp;
   30:     stemp = sdot = 0.0;
   31:     if(n <= 0) return sdot;
   32:     if(incx == 1 && incy == 1) goto S20;
   33:     ix = iy = 1;
   34:     if(incx < 0) ix = (-n+1)*incx+1;
   35:     if(incy < 0) iy = (-n+1)*incy+1;
   36:     for(i=1; i<=n; i++) {
   37:         stemp += (*(sx+ix-1)**(sy+iy-1));
   38:         ix += incx;
   39:         iy += incy;
   40:     }
   41:     sdot = stemp;
   42:     return sdot;
   43: S20:
   44:     m = n % 5L;
   45:     if(m == 0) goto S40;
   46:     for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
   47:     if(n < 5) goto S60;
   48: S40:
   49:     mp1 = m+1;
   50:     for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
   51:       +1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
   52: S60:
   53:     sdot = stemp;
   54:     return sdot;
   55: }
   56: void spofa(float *a,long lda,long n,long *info)
   57: /*
   58:      SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
   59:      SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
   60:      DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
   61:      (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
   62:      ON ENTRY
   63:         A       REAL(LDA, N)
   64:                 THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
   65:                 DIAGONAL AND UPPER TRIANGLE ARE USED.
   66:         LDA     INTEGER
   67:                 THE LEADING DIMENSION OF THE ARRAY  A .
   68:         N       INTEGER
   69:                 THE ORDER OF THE MATRIX  A .
   70:      ON RETURN
   71:         A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
   72:                 WHERE  TRANS(R)  IS THE TRANSPOSE.
   73:                 THE STRICT LOWER TRIANGLE IS UNALTERED.
   74:                 IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
   75:         INFO    INTEGER
   76:                 = 0  FOR NORMAL RETURN.
   77:                 = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
   78:                      OF ORDER  K  IS NOT POSITIVE DEFINITE.
   79:      LINPACK.  THIS VERSION DATED 08/14/78 .
   80:      CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
   81:      SUBROUTINES AND FUNCTIONS
   82:      BLAS SDOT
   83:      FORTRAN SQRT
   84:      INTERNAL VARIABLES
   85: */
   86: {
   87: extern float sdot(long n,float *sx,long incx,float *sy,long incy);
   88: static long j,jm1,k;
   89: static float t,s;
   90: /*
   91:      BEGIN BLOCK WITH ...EXITS TO 40
   92: */
   93:     for(j=1; j<=n; j++) {
   94:         *info = j;
   95:         s = 0.0;
   96:         jm1 = j-1;
   97:         if(jm1 < 1) goto S20;
   98:         for(k=0; k<jm1; k++) {
   99:             t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
  100:             t /=  *(a+k+k*lda);
  101:             *(a+k+(j-1)*lda) = t;
  102:             s += (t*t);
  103:         }
  104: S20:
  105:         s = *(a+j-1+(j-1)*lda)-s;
  106: /*
  107:      ......EXIT
  108: */
  109:         if(s <= 0.0) goto S40;
  110:         *(a+j-1+(j-1)*lda) = sqrt(s);
  111:     }
  112:     *info = 0;
  113: S40:
  114:     return;
  115: }