[Bio] / Numeric / Src / linpack.c Repository:
ViewVC logotype

Annotation of /Numeric/Src/linpack.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (view) (download) (as text)

1 : efrank 1.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 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3