[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / d / linsq.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:21  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE LINSQ(K,N,M,F,X,Y,W,DA,H,COV,QPRINT,QMIN)
      DIMENSION X(M,N),F(K),Y(N),W(N),H(K,K),DA(K)
      DIMENSION INDEX(100)
      DOUBLE PRECISION H,DETERM
      DOUBLE PRECISION DWI,DYI,DFK1
      KPLUS=K+1
      APRINT=QPRINT
      DO 11 I=1,K
      DO 11 J=1,KPLUS
   11 H(I,J) = 0.D0
      IFLAG=1
      DO 2 I=1,N
      CALL FCN(K,M,F,X(1,I),IFLAG)
      IFLAG=4
      DWI = DBLE(W(I))
      DYI = DBLE(Y(I))
      DO 2 K1=1,K
      DFK1 = DBLE(F(K1))
      H(K1,K+1) = H(K1,K+1) + DWI*DFK1*DYI
      DO 2 L=1,K1
    2 H(K1,L) = H(K1,L) + DWI*DFK1*DBLE(F(L))
      DO 3 J=1,K
      DO 3 I=1,J
    3 H(I,J)=H(J,I)
      IF(COV)5,4,5
    4 CALL D508R2(H,K,K,K+1,1,INDEX,NERROR,DETERM)
      MM=1
      QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
      DO 7 I=1,K
    7 F(I)=H(I,1)
      IF(APRINT.NE.0.) GO TO 10
      RETURN
    5 CALL D508R1(H,K,K,K+1,1,INDEX,NERROR,DETERM)
      MM=K+1
      QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
      DO 6 I=1,K
    6 F(I)=H(I,K+1)
      DO 8 I=1,K
      HDA=H(I,I)
    8 DA(I)=SQRT(HDA)
      IF(APRINT.EQ.0.) RETURN
C
C     PRINTS ERRORS IN THE COEFFICIENTS
C
      WRITE(6,100)
  100 FORMAT(45X,'ERRORS IN THE COEFFICIENTS')
      WRITE(6,101)(DA(I),I=1,K)
  101 FORMAT(5X,5G16.6)
C
C     PRINTS LOWER DIAGONAL COVARIANCE MATRIX
C
      WRITE(6,102)
  102 FORMAT(///40X,'VARIANCE-COVARIANCE MATRIX'///)
      DO 9 I=1,K
      DO 12 J=1,K
   12 F(J)=H(I,J)
    9 WRITE(6,103)(F(JJ),JJ=1,I)
  103 FORMAT(10X,5G20.6//)
      DO 13 I=1,K
   13 F(I)=H(I,K+1)
   10 CONTINUE
C
C     PRINTS COEFFICIENTS AND SUM OF SQUARES
C
      WRITE(6,104)QMIN
  104 FORMAT(///30X,'SUM OF SQUARES=',G20.6///)
      WRITE(6,105)
  105 FORMAT(45X,'COEFFICIENTS'///)
      WRITE(6,103)(F(I),I=1,K)
      RETURN
      END
      FUNCTION QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
C     COMPUTES THE SUM OF SQUARES
      DIMENSION F(K),X(M,N),DA(K),H(K,K),Y(N),W(N)
      DOUBLE PRECISION H
      IFLAG = 4
      Q=0.
      DO 3 I=1,N
      CALL FCN(K,M,F,X(1,I),IFLAG)
      HSUM=0.
      DO 2 I1=1,K
    2 HSUM=H(I1,MM)*F(I1)+HSUM
    3 Q=(HSUM-Y(I))*(HSUM-Y(I))*W(I)+Q
      QLINSQ = Q
      RETURN
      END
      SUBROUTINE  D508R1 (A,IDIM1,N1,IDIM2,N2,INDEX,NERROR,DETERM)
C
C        MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS
      DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
      DIMENSION A(IDIM1),INDEX(IDIM1)
      DETER=1.0D0
      N=N1
      IEMAT=N+N2
      IDIM=IDIM1
      NMIN1=N-1
C        THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
C        ARRAY A.
      IPIVC=1-IDIM
C        MAIN LOOP TO INVERT THE MATRIX
      DO 11 MAIN=1,N
      PIVOT=0.0D0
      IPIVC=IPIVC+IDIM
C        SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
      IPIVC1=IPIVC+MAIN-1
      IPIVC2=IPIVC +NMIN1
      DO 2 I1=IPIVC1,IPIVC2
      IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
    1 PIVOT=A(I1)
      LPIV=I1
    2 CONTINUE
C        IS PIVOT DIFFERENT FROM ZERO
      IF(PIVOT) 3,15,3
C        GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
    3 ICOL=LPIV-IPIVC+1
      INDEX(MAIN)=ICOL
      IF(ICOL-MAIN) 6,6,4
C        COMPLEMENT THE DETERMINANT
    4 DETER=-DETER
C        POINTER TO LINE PIVOT FOUND
      ICOL=ICOL-IDIM
C        POINTER TO EXACT PIVOT LINE
      I3=MAIN-IDIM
      DO 5 I=1,IEMAT
      ICOL=ICOL+IDIM
      I3=I3+IDIM
      SWAP=A(I3)
      A(I3)=A(ICOL)
    5 A(ICOL)=SWAP
C        COMPUTE DETERMINANT
    6 DETER=DETER*PIVOT
      PIVOT=1./PIVOT
C        TRANSFORM PIVOT COLUMN
      I3=IPIVC+NMIN1
      DO 7 I=IPIVC,I3
    7 A(I)=-A(I)*PIVOT
      A(IPIVC1)=PIVOT
C        PIVOT ELEMENT TRANSFORMED
C
C        NOW CONVERT REST OF THE MATRIX
      I1=MAIN-IDIM
C        POINTER TO PIVOT LINE ELEMENTS
      ICOL=1-IDIM
C        GENERAL COLUMN POINTER
      DO 10 I=1,IEMAT
      ICOL=ICOL+IDIM
      I1=I1+IDIM
C        POINTERS MOVED
      IF(I-MAIN) 8,10,8
C        PIVOT COLUMN EXCLUDED
    8 JCOL=ICOL+NMIN1
      SWAP=A(I1)
      I3=IPIVC-1
      DO 9 I2=ICOL,JCOL
      I3=I3+1
    9 A(I2)=A(I2)+SWAP*A(I3)
      A(I1)=SWAP*PIVOT
   10 CONTINUE
   11 CONTINUE
C        NOW REARRANGE THE MATRIX TO GET RIGHT INVERS
      DO 14 I1=1,N
      MAIN=N+1-I1
      LPIV=INDEX(MAIN)
      IF(LPIV-MAIN) 12,14,12
   12 ICOL=(LPIV-1)*IDIM+1
      JCOL=ICOL+NMIN1
      IPIVC=(MAIN-1)*IDIM+1-ICOL
      DO 13 I2=ICOL,JCOL
      I3=I2+IPIVC
      SWAP=A(I2)
      A(I2)=A(I3)
   13 A(I3)=SWAP
   14 CONTINUE
      DETERM=DETER
      NERROR=0
      RETURN
   15 NERROR=MAIN
      DETERM=DETER
      RETURN
      END
      SUBROUTINE D508R2 (A,DIM1,N1,DIM2,N2,INDEX,NERROR,DETERM)
      INTEGER DIM1,DIM,PIVCOL,PIVCO1,TOPX,ENDX,TOPCOL,ENDCOL,EMAT
      DIMENSION A(DIM1),INDEX(DIM1)
      DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
      DIM=DIM1
      DETER=1.0D0
      N=N1
      EMAT=N+N2
      NMIN1=N-1
      PIVCOL=-DIM
C     MAIN LOOP TO CREATE TRIANGULAR
      DO 10 MAIN=1,N
      PIVOT=0.0D0
      PIVCOL=PIVCOL+DIM+1
      PIVCO1=PIVCOL+N-MAIN
C     SEARCH PIVOT
      DO 2 I1=PIVCOL,PIVCO1
      IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
    1 PIVOT=A(I1)
      LPIV=I1
    2 CONTINUE
C     IS PIVOT DIFFERENT FROM ZERO
      IF(PIVOT) 3,15,3
C     IS IT NECESSARY TO BRING PIVOT TO DIAGONAL
    3 IF(LPIV-PIVCOL) 4,6,4
    4 DETER=-DETER
      LPIV=LPIV-DIM
      I1=PIVCOL-DIM
      DO 5 I2=MAIN,EMAT
      LPIV=LPIV+DIM
      I1=I1+DIM
      SWAP=A(I1)
      A(I1)=A(LPIV)
    5 A(LPIV)=SWAP
    6 DETER=DETER*PIVOT
      IF (MAIN .EQ. N)  GO TO 10
      PIVOT=1./PIVOT
C     MODIFY PIVOT COLUMN
      I1=PIVCOL+1
      DO 7 I2=I1,PIVCO1
    7 A(I2)=A(I2)*PIVOT
C     CONVERT THE SUBMATRIX AND RIGHT SIDES
      I3=PIVCOL
      IROW=MAIN+1
      DO 9 I1=IROW,N
      I3=I3+1
      I4=PIVCOL
      I5=I3
      DO 8 I2=IROW,EMAT
      I4=I4+DIM
      I5=I5+DIM
    8 A(I5)=A(I5)-A(I4)*A(I3)
    9 CONTINUE
   10 CONTINUE
      DETERM=DETER
      NERROR=0
C     COMPUTE THE SOLUTIONS
      NO=N+1
      TOPX=NMIN1*DIM+1
      DO 13 I=NO,EMAT
      TOPX=TOPX+DIM
      ENDX=TOPX+N
      TOPCOL=N*DIM+1
      ENDCOL=TOPCOL+NMIN1
      DO 12 I1=1,NMIN1
      ENDX=ENDX-1
      TOPCOL=TOPCOL-DIM
      ENDCOL=ENDCOL-DIM-1
      A(ENDX)=A(ENDX)/A(ENDCOL+1)
      SWAP=A(ENDX)
      I3=TOPX-1
      DO 11 I2=TOPCOL,ENDCOL
      I3=I3+1
   11 A(I3)=A(I3)-A(I2)*SWAP
   12 CONTINUE
      A(TOPX)=A(TOPX)/A(1)
   13 CONTINUE
C     LEFTADJUST THE SOLUTIONS
      I=-DIM
      TOPX=NMIN1*DIM+1
      ENDX=TOPX+NMIN1
      DO 14 I1=NO,EMAT
      TOPX=TOPX+DIM
      ENDX=ENDX+DIM
      I=I+DIM
      I3=I
      DO 14 I2=TOPX,ENDX
      I3=I3+1
   14 A(I3)=A(I2)
      RETURN
C     ERROR EXIT
   15 NERROR=-1
      DETERM=DETER
      WRITE(6,100)MAIN,MAIN
      RETURN
  100 FORMAT(' LINEQ1 ..... THE ',I10,'. COLUMN OF THE MATRIX CONTAINS'
     1 ,' ONLY ZEROS AT THE',I10,'. ELIMINATIONSTEP')
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:21 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	SUBROUTINE LINSQ(K,N,M,F,X,Y,W,DA,H,COV,QPRINT,QMIN)
	11	DIMENSION X(M,N),F(K),Y(N),W(N),H(K,K),DA(K)
	12	DIMENSION INDEX(100)
	13	DOUBLE PRECISION H,DETERM
	14	DOUBLE PRECISION DWI,DYI,DFK1
	15	KPLUS=K+1
	16	APRINT=QPRINT
	17	DO 11 I=1,K
	18	DO 11 J=1,KPLUS
	19	11 H(I,J) = 0.D0
	20	IFLAG=1
	21	DO 2 I=1,N
	22	CALL FCN(K,M,F,X(1,I),IFLAG)
	23	IFLAG=4
	24	DWI = DBLE(W(I))
	25	DYI = DBLE(Y(I))
	26	DO 2 K1=1,K
	27	DFK1 = DBLE(F(K1))
	28	H(K1,K+1) = H(K1,K+1) + DWIDFK1DYI
	29	DO 2 L=1,K1
	30	2 H(K1,L) = H(K1,L) + DWIDFK1DBLE(F(L))
	31	DO 3 J=1,K
	32	DO 3 I=1,J
	33	3 H(I,J)=H(J,I)
	34	IF(COV)5,4,5
	35	4 CALL D508R2(H,K,K,K+1,1,INDEX,NERROR,DETERM)
	36	MM=1
	37	QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
	38	DO 7 I=1,K
	39	7 F(I)=H(I,1)
	40	IF(APRINT.NE.0.) GO TO 10
	41	RETURN
	42	5 CALL D508R1(H,K,K,K+1,1,INDEX,NERROR,DETERM)
	43	MM=K+1
	44	QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
	45	DO 6 I=1,K
	46	6 F(I)=H(I,K+1)
	47	DO 8 I=1,K
	48	HDA=H(I,I)
	49	8 DA(I)=SQRT(HDA)
	50	IF(APRINT.EQ.0.) RETURN
	51	C
	52	C PRINTS ERRORS IN THE COEFFICIENTS
	53	C
	54	WRITE(6,100)
	55	100 FORMAT(45X,'ERRORS IN THE COEFFICIENTS')
	56	WRITE(6,101)(DA(I),I=1,K)
	57	101 FORMAT(5X,5G16.6)
	58	C
	59	C PRINTS LOWER DIAGONAL COVARIANCE MATRIX
	60	C
	61	WRITE(6,102)
	62	102 FORMAT(///40X,'VARIANCE-COVARIANCE MATRIX'///)
	63	DO 9 I=1,K
	64	DO 12 J=1,K
65	12 F(J)=H(I,J)
66	9 WRITE(6,103)(F(JJ),JJ=1,I)
67	103 FORMAT(10X,5G20.6//)
68	DO 13 I=1,K
69	13 F(I)=H(I,K+1)
70	10 CONTINUE
71	C
72	C PRINTS COEFFICIENTS AND SUM OF SQUARES
73	C
74	WRITE(6,104)QMIN
75	104 FORMAT(///30X,'SUM OF SQUARES=',G20.6///)
76	WRITE(6,105)
77	105 FORMAT(45X,'COEFFICIENTS'///)
78	WRITE(6,103)(F(I),I=1,K)
79	RETURN
80	END
81	FUNCTION QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
82	C COMPUTES THE SUM OF SQUARES
83	DIMENSION F(K),X(M,N),DA(K),H(K,K),Y(N),W(N)
84	DOUBLE PRECISION H
85	IFLAG = 4
86	Q=0.
87	DO 3 I=1,N
88	CALL FCN(K,M,F,X(1,I),IFLAG)
89	HSUM=0.
90	DO 2 I1=1,K
91	2 HSUM=H(I1,MM)*F(I1)+HSUM
92	3 Q=(HSUM-Y(I))(HSUM-Y(I))W(I)+Q
93	QLINSQ = Q
94	RETURN
95	END
96	SUBROUTINE D508R1 (A,IDIM1,N1,IDIM2,N2,INDEX,NERROR,DETERM)
97	C
98	C MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS
99	DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
100	DIMENSION A(IDIM1),INDEX(IDIM1)
101	DETER=1.0D0
102	N=N1
103	IEMAT=N+N2
104	IDIM=IDIM1
105	NMIN1=N-1
106	C THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
107	C ARRAY A.
108	IPIVC=1-IDIM
109	C MAIN LOOP TO INVERT THE MATRIX
110	DO 11 MAIN=1,N
111	PIVOT=0.0D0
112	IPIVC=IPIVC+IDIM
113	C SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
114	IPIVC1=IPIVC+MAIN-1
115	IPIVC2=IPIVC +NMIN1
116	DO 2 I1=IPIVC1,IPIVC2
117	IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
118	1 PIVOT=A(I1)
119	LPIV=I1
120	2 CONTINUE
121	C IS PIVOT DIFFERENT FROM ZERO
122	IF(PIVOT) 3,15,3
123	C GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
124	3 ICOL=LPIV-IPIVC+1
125	INDEX(MAIN)=ICOL
126	IF(ICOL-MAIN) 6,6,4
127	C COMPLEMENT THE DETERMINANT
128	4 DETER=-DETER
129	C POINTER TO LINE PIVOT FOUND
130	ICOL=ICOL-IDIM
131	C POINTER TO EXACT PIVOT LINE
132	I3=MAIN-IDIM
133	DO 5 I=1,IEMAT
134	ICOL=ICOL+IDIM
135	I3=I3+IDIM
136	SWAP=A(I3)
137	A(I3)=A(ICOL)
138	5 A(ICOL)=SWAP
139	C COMPUTE DETERMINANT
140	6 DETER=DETER*PIVOT
141	PIVOT=1./PIVOT
142	C TRANSFORM PIVOT COLUMN
143	I3=IPIVC+NMIN1
144	DO 7 I=IPIVC,I3
145	7 A(I)=-A(I)*PIVOT
146	A(IPIVC1)=PIVOT
147	C PIVOT ELEMENT TRANSFORMED
148	C
149	C NOW CONVERT REST OF THE MATRIX
150	I1=MAIN-IDIM
151	C POINTER TO PIVOT LINE ELEMENTS
152	ICOL=1-IDIM
153	C GENERAL COLUMN POINTER
154	DO 10 I=1,IEMAT
155	ICOL=ICOL+IDIM
156	I1=I1+IDIM
157	C POINTERS MOVED
158	IF(I-MAIN) 8,10,8
159	C PIVOT COLUMN EXCLUDED
160	8 JCOL=ICOL+NMIN1
161	SWAP=A(I1)
162	I3=IPIVC-1
163	DO 9 I2=ICOL,JCOL
164	I3=I3+1
165	9 A(I2)=A(I2)+SWAP*A(I3)
166	A(I1)=SWAP*PIVOT
167	10 CONTINUE
168	11 CONTINUE
169	C NOW REARRANGE THE MATRIX TO GET RIGHT INVERS
170	DO 14 I1=1,N
171	MAIN=N+1-I1
172	LPIV=INDEX(MAIN)
173	IF(LPIV-MAIN) 12,14,12
174	12 ICOL=(LPIV-1)*IDIM+1
175	JCOL=ICOL+NMIN1
176	IPIVC=(MAIN-1)*IDIM+1-ICOL
177	DO 13 I2=ICOL,JCOL
178	I3=I2+IPIVC
179	SWAP=A(I2)
180	A(I2)=A(I3)
181	13 A(I3)=SWAP
182	14 CONTINUE
183	DETERM=DETER
184	NERROR=0
185	RETURN
186	15 NERROR=MAIN
187	DETERM=DETER
188	RETURN
189	END
190	SUBROUTINE D508R2 (A,DIM1,N1,DIM2,N2,INDEX,NERROR,DETERM)
191	INTEGER DIM1,DIM,PIVCOL,PIVCO1,TOPX,ENDX,TOPCOL,ENDCOL,EMAT
192	DIMENSION A(DIM1),INDEX(DIM1)
193	DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
194	DIM=DIM1
195	DETER=1.0D0
196	N=N1
197	EMAT=N+N2
198	NMIN1=N-1
199	PIVCOL=-DIM
200	C MAIN LOOP TO CREATE TRIANGULAR
201	DO 10 MAIN=1,N
202	PIVOT=0.0D0
203	PIVCOL=PIVCOL+DIM+1
204	PIVCO1=PIVCOL+N-MAIN
205	C SEARCH PIVOT
206	DO 2 I1=PIVCOL,PIVCO1
207	IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
208	1 PIVOT=A(I1)
209	LPIV=I1
210	2 CONTINUE
211	C IS PIVOT DIFFERENT FROM ZERO
212	IF(PIVOT) 3,15,3
213	C IS IT NECESSARY TO BRING PIVOT TO DIAGONAL
214	3 IF(LPIV-PIVCOL) 4,6,4
215	4 DETER=-DETER
216	LPIV=LPIV-DIM
217	I1=PIVCOL-DIM
218	DO 5 I2=MAIN,EMAT
219	LPIV=LPIV+DIM
220	I1=I1+DIM
221	SWAP=A(I1)
222	A(I1)=A(LPIV)
223	5 A(LPIV)=SWAP
224	6 DETER=DETER*PIVOT
225	IF (MAIN .EQ. N) GO TO 10
226	PIVOT=1./PIVOT
227	C MODIFY PIVOT COLUMN
228	I1=PIVCOL+1
229	DO 7 I2=I1,PIVCO1
230	7 A(I2)=A(I2)*PIVOT
231	C CONVERT THE SUBMATRIX AND RIGHT SIDES
232	I3=PIVCOL
233	IROW=MAIN+1
234	DO 9 I1=IROW,N
235	I3=I3+1
236	I4=PIVCOL
237	I5=I3
238	DO 8 I2=IROW,EMAT
239	I4=I4+DIM
240	I5=I5+DIM
241	8 A(I5)=A(I5)-A(I4)*A(I3)
242	9 CONTINUE
243	10 CONTINUE
244	DETERM=DETER
245	NERROR=0
246	C COMPUTE THE SOLUTIONS
247	NO=N+1
248	TOPX=NMIN1*DIM+1
249	DO 13 I=NO,EMAT
250	TOPX=TOPX+DIM
251	ENDX=TOPX+N
252	TOPCOL=N*DIM+1
253	ENDCOL=TOPCOL+NMIN1
254	DO 12 I1=1,NMIN1
255	ENDX=ENDX-1
256	TOPCOL=TOPCOL-DIM
257	ENDCOL=ENDCOL-DIM-1
258	A(ENDX)=A(ENDX)/A(ENDCOL+1)
259	SWAP=A(ENDX)
260	I3=TOPX-1
261	DO 11 I2=TOPCOL,ENDCOL
262	I3=I3+1
263	11 A(I3)=A(I3)-A(I2)*SWAP
264	12 CONTINUE
265	A(TOPX)=A(TOPX)/A(1)
266	13 CONTINUE
267	C LEFTADJUST THE SOLUTIONS
268	I=-DIM
269	TOPX=NMIN1*DIM+1
270	ENDX=TOPX+NMIN1
271	DO 14 I1=NO,EMAT
272	TOPX=TOPX+DIM
273	ENDX=ENDX+DIM
274	I=I+DIM
275	I3=I
276	DO 14 I2=TOPX,ENDX
277	I3=I3+1
278	14 A(I3)=A(I2)
279	RETURN
280	C ERROR EXIT
281	15 NERROR=-1
282	DETERM=DETER
283	WRITE(6,100)MAIN,MAIN
284	RETURN
285	100 FORMAT(' LINEQ1 ..... THE ',I10,'. COLUMN OF THE MATRIX CONTAINS'
286	1 ,' ONLY ZEROS AT THE',I10,'. ELIMINATIONSTEP')
287	END