| line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
|
1
|
|
|
|
|
|
|
#include |
|
2
|
3
|
|
|
|
|
|
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
|
3
|
|
|
|
|
|
stemp = sdot = 0.0; |
|
7
|
3
|
50
|
|
|
|
|
if(n <= 0) return sdot; |
|
8
|
0
|
0
|
|
|
|
|
if(incx == 1 && incy == 1) goto S20; |
|
|
|
0
|
|
|
|
|
|
|
9
|
0
|
|
|
|
|
|
ix = iy = 1; |
|
10
|
0
|
0
|
|
|
|
|
if(incx < 0) ix = (-n+1)*incx+1; |
|
11
|
0
|
0
|
|
|
|
|
if(incy < 0) iy = (-n+1)*incy+1; |
|
12
|
0
|
0
|
|
|
|
|
for(i=1; i<=n; i++) { |
|
13
|
0
|
|
|
|
|
|
stemp += (*(sx+ix-1)**(sy+iy-1)); |
|
14
|
0
|
|
|
|
|
|
ix += incx; |
|
15
|
0
|
|
|
|
|
|
iy += incy; |
|
16
|
|
|
|
|
|
|
} |
|
17
|
0
|
|
|
|
|
|
sdot = stemp; |
|
18
|
0
|
|
|
|
|
|
return sdot; |
|
19
|
0
|
|
|
|
|
|
S20: |
|
20
|
0
|
|
|
|
|
|
m = n % 5L; |
|
21
|
0
|
0
|
|
|
|
|
if(m == 0) goto S40; |
|
22
|
0
|
0
|
|
|
|
|
for(i=0; i
|
|
23
|
0
|
0
|
|
|
|
|
if(n < 5) goto S60; |
|
24
|
0
|
|
|
|
|
|
S40: |
|
25
|
0
|
|
|
|
|
|
mp1 = m+1; |
|
26
|
0
|
0
|
|
|
|
|
for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i |
|
27
|
0
|
|
|
|
|
|
+1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3)); |
|
28
|
0
|
|
|
|
|
|
S60: |
|
29
|
0
|
|
|
|
|
|
sdot = stemp; |
|
30
|
0
|
|
|
|
|
|
return sdot; |
|
31
|
|
|
|
|
|
|
} |
|
32
|
3
|
|
|
|
|
|
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
|
9
|
100
|
|
|
|
|
for(j=1; j<=n; j++) { |
|
70
|
6
|
|
|
|
|
|
*info = j; |
|
71
|
6
|
|
|
|
|
|
s = 0.0; |
|
72
|
6
|
|
|
|
|
|
jm1 = j-1; |
|
73
|
6
|
100
|
|
|
|
|
if(jm1 < 1) goto S20; |
|
74
|
6
|
100
|
|
|
|
|
for(k=0; k
|
|
75
|
3
|
|
|
|
|
|
t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L); |
|
76
|
3
|
|
|
|
|
|
t /= *(a+k+k*lda); |
|
77
|
3
|
|
|
|
|
|
*(a+k+(j-1)*lda) = t; |
|
78
|
3
|
|
|
|
|
|
s += (t*t); |
|
79
|
|
|
|
|
|
|
} |
|
80
|
3
|
|
|
|
|
|
S20: |
|
81
|
6
|
|
|
|
|
|
s = *(a+j-1+(j-1)*lda)-s; |
|
82
|
|
|
|
|
|
|
/* |
|
83
|
|
|
|
|
|
|
......EXIT |
|
84
|
|
|
|
|
|
|
*/ |
|
85
|
6
|
50
|
|
|
|
|
if(s <= 0.0) goto S40; |
|
86
|
6
|
|
|
|
|
|
*(a+j-1+(j-1)*lda) = sqrt(s); |
|
87
|
|
|
|
|
|
|
} |
|
88
|
3
|
|
|
|
|
|
*info = 0; |
|
89
|
3
|
|
|
|
|
|
S40: |
|
90
|
3
|
|
|
|
|
|
return; |
|
91
|
|
|
|
|
|
|
} |