File:
[LON-CAPA] /
capa /
capa51 /
pProj /
linpack.c
Revision
1.4:
download - view:
text,
annotated -
select for diffs
Mon Aug 7 20:47:29 2000 UTC (24 years, 1 month ago) by
albertel
Branches:
MAIN
CVS tags:
version_2_9_X,
version_2_9_99_0,
version_2_9_1,
version_2_9_0,
version_2_8_X,
version_2_8_99_1,
version_2_8_99_0,
version_2_8_2,
version_2_8_1,
version_2_8_0,
version_2_7_X,
version_2_7_99_1,
version_2_7_99_0,
version_2_7_1,
version_2_7_0,
version_2_6_X,
version_2_6_99_1,
version_2_6_99_0,
version_2_6_3,
version_2_6_2,
version_2_6_1,
version_2_6_0,
version_2_5_X,
version_2_5_99_1,
version_2_5_99_0,
version_2_5_2,
version_2_5_1,
version_2_5_0,
version_2_4_X,
version_2_4_99_0,
version_2_4_2,
version_2_4_1,
version_2_4_0,
version_2_3_X,
version_2_3_99_0,
version_2_3_2,
version_2_3_1,
version_2_3_0,
version_2_2_X,
version_2_2_99_1,
version_2_2_99_0,
version_2_2_2,
version_2_2_1,
version_2_2_0,
version_2_1_X,
version_2_1_99_3,
version_2_1_99_2,
version_2_1_99_1,
version_2_1_99_0,
version_2_1_3,
version_2_1_2,
version_2_1_1,
version_2_1_0,
version_2_12_X,
version_2_11_X,
version_2_11_5_msu,
version_2_11_5,
version_2_11_4_uiuc,
version_2_11_4_msu,
version_2_11_4,
version_2_11_3_uiuc,
version_2_11_3_msu,
version_2_11_3,
version_2_11_2_uiuc,
version_2_11_2_msu,
version_2_11_2_educog,
version_2_11_2,
version_2_11_1,
version_2_11_0_RC3,
version_2_11_0_RC2,
version_2_11_0_RC1,
version_2_11_0,
version_2_10_X,
version_2_10_1,
version_2_10_0_RC2,
version_2_10_0_RC1,
version_2_10_0,
version_2_0_X,
version_2_0_99_1,
version_2_0_2,
version_2_0_1,
version_2_0_0,
version_1_99_3,
version_1_99_2,
version_1_99_1_tmcc,
version_1_99_1,
version_1_99_0_tmcc,
version_1_99_0,
version_1_3_X,
version_1_3_3,
version_1_3_2,
version_1_3_1,
version_1_3_0,
version_1_2_X,
version_1_2_99_1,
version_1_2_99_0,
version_1_2_1,
version_1_2_0,
version_1_1_X,
version_1_1_99_5,
version_1_1_99_4,
version_1_1_99_3,
version_1_1_99_2,
version_1_1_99_1,
version_1_1_99_0,
version_1_1_3,
version_1_1_2,
version_1_1_1,
version_1_1_0,
version_1_0_99_3,
version_1_0_99_2,
version_1_0_99_1,
version_1_0_99,
version_1_0_3,
version_1_0_2,
version_1_0_1,
version_1_0_0,
version_0_99_5,
version_0_99_4,
version_0_99_3,
version_0_99_2,
version_0_99_1,
version_0_99_0,
version_0_6_2,
version_0_6,
version_0_5_1,
version_0_5,
version_0_4,
stable_2002_spring,
stable_2002_july,
stable_2002_april,
stable_2001_fall,
release_5-1-3,
loncapaMITrelate_1,
language_hyphenation_merge,
language_hyphenation,
conference_2003,
bz6209-base,
bz6209,
STABLE,
HEAD,
GCI_3,
GCI_2,
GCI_1,
CAPA_5-1-6,
CAPA_5-1-5,
CAPA_5-1-4_RC1,
BZ4492-merge,
BZ4492-feature_horizontal_radioresponse,
BZ4492-feature_Support_horizontal_radioresponse,
BZ4492-Support_horizontal_radioresponse
- fixed license notices the reference the GNU GPL rather than the GNU LGPL
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 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: General Public License for more details.
13:
14: You should have received a copy of the GNU 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>