1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#include <math.h>
double sdot(long n,double *sx,long incx,double *sy,long incy)
{
static long i,ix,iy,m,mp1;
static double sdot,stemp;
stemp = sdot = 0.0;
if(n <= 0) return sdot;
if(incx == 1 && incy == 1) goto S20;
ix = iy = 1;
if(incx < 0) ix = (-n+1)*incx+1;
if(incy < 0) iy = (-n+1)*incy+1;
for(i=1; i<=n; i++) {
stemp += (*(sx+ix-1)**(sy+iy-1));
ix += incx;
iy += incy;
}
sdot = stemp;
return sdot;
S20:
m = n % 5L;
if(m == 0) goto S40;
for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
if(n < 5) goto S60;
S40:
mp1 = m+1;
for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
+1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
S60:
sdot = stemp;
return sdot;
}
void spofa(double *a,long lda,long n,long *info)
/*
SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
(TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
ON ENTRY
A REAL(LDA, N)
THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE
DIAGONAL AND UPPER TRIANGLE ARE USED.
LDA INTEGER
THE LEADING DIMENSION OF THE ARRAY A .
N INTEGER
THE ORDER OF THE MATRIX A .
ON RETURN
A AN UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R
WHERE TRANS(R) IS THE TRANSPOSE.
THE STRICT LOWER TRIANGLE IS UNALTERED.
IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
INFO INTEGER
= 0 FOR NORMAL RETURN.
= K SIGNALS AN ERROR CONDITION. THE LEADING MINOR
OF ORDER K IS NOT POSITIVE DEFINITE.
LINPACK. THIS VERSION DATED 08/14/78 .
CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
SUBROUTINES AND FUNCTIONS
BLAS SDOT
FORTRAN SQRT
INTERNAL VARIABLES
*/
{
extern double sdot(long n,double *sx,long incx,double *sy,long incy);
static long j,jm1,k;
static double t,s;
/*
BEGIN BLOCK WITH ...EXITS TO 40
*/
for(j=1; j<=n; j++) {
*info = j;
s = 0.0;
jm1 = j-1;
if(jm1 < 1) goto S20;
for(k=0; k<jm1; k++) {
t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
t /= *(a+k+k*lda);
*(a+k+(j-1)*lda) = t;
s += (t*t);
}
S20:
s = *(a+j-1+(j-1)*lda)-s;
/*
......EXIT
*/
if(s <= 0.0) goto S40;
*(a+j-1+(j-1)*lda) = sqrt(s);
}
*info = 0;
S40:
return;
}
|
