modules/Math-Random/linpack.c

   1 #include <math.h>
   2 double sdot(long n,double *sx,long incx,double *sy,long incy)
   3 {
   4 static long i,ix,iy,m,mp1;
   5 static double 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(double *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 double sdot(long n,double *sx,long incx,double *sy,long incy);
  64 static long j,jm1,k;
  65 static double 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 }