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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
*DECK LIBOMG
SUBROUTINE LIBOMG(MX,IX,X,MY,IY,Y,DCM,OMEG)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the correlated weight matrix preserving a matrix of moments.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* MX number of base points in the first variable.
* IX order of the first moment of the first variable. We must
* have 1-MX <= IX <= 0 (order 0 moment is always preserved).
* X base points in the first variable.
* MY number of base points in the second variable.
* IY order of the first moment of the second variable. We must
* have 1-MY <= IY <= 0 (order 0 moment is always preserved).
* Y base points in the second variable.
* DCM co-moments.
*
*Parameters: output
* OMEG correlated weight matrix.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
REAL X(MX),Y(MY),OMEG(MX,MY)
DOUBLE PRECISION DCM(MX,MY)
*----
* LOCAL VARIABLES
*----
PARAMETER (MAXNOR=20)
DOUBLE PRECISION DD,DAUX,WORK,DDA(0:MAXNOR),PROD1(MAXNOR,MAXNOR),
1 PROD2(MAXNOR,MAXNOR)
*
IF(MX.GT.MAXNOR) CALL XABORT('LIBOMG: STORAGE OVERFLOW(1).')
IF(MY.GT.MAXNOR) CALL XABORT('LIBOMG: STORAGE OVERFLOW(2).')
IF((1-MX.GT.IX).OR.(IX.GT.0)) CALL XABORT('LIBOMG: INCONSISTENT '
1 //'VALUE OF IX.')
IF((1-MY.GT.IY).OR.(IY.GT.0)) CALL XABORT('LIBOMG: INCONSISTENT '
1 //'VALUE OF IY.')
*
DO 15 I=1,MX
DO 10 J=1,MY
PROD1(I,J)=0.0D0
10 CONTINUE
15 CONTINUE
DO 52 I=1,MX
DAUX=DBLE(X(I))
DDA(0)=1.0D0
J0=0
DO 30 J=1,MX
IF(J.EQ.I) GO TO 30
J0=J0+1
DDA(J0)=DDA(J0-1)
DO 20 K=1,J0-1
DDA(J0-K)=DDA(J0-K-1)-DDA(J0-K)*DBLE(X(J))
20 CONTINUE
DDA(0)=-DDA(0)*DBLE(X(J))
30 CONTINUE
DD=1.0D0
DO 40 J=1,MX
IF(J.NE.I) DD=DD*(DBLE(X(J))-DAUX)
40 CONTINUE
WORK=((-1.0D0)**(MX-1))*DAUX**(-IX)/DD
DO 51 J=1,MY
DO 50 K=1,MX
PROD1(I,J)=PROD1(I,J)+WORK*DDA(K-1)*DCM(K,J)
50 CONTINUE
51 CONTINUE
52 CONTINUE
*
DO 65 I=1,MX
DO 60 J=1,MY
PROD2(I,J)=0.0D0
60 CONTINUE
65 CONTINUE
DO 102 I=1,MY
DAUX=DBLE(Y(I))
DDA(0)=1.0D0
J0=0
DO 80 J=1,MY
IF(J.EQ.I) GO TO 80
J0=J0+1
DDA(J0)=DDA(J0-1)
DO 70 K=1,J0-1
DDA(J0-K)=DDA(J0-K-1)-DDA(J0-K)*DBLE(Y(J))
70 CONTINUE
DDA(0)=-DDA(0)*DBLE(Y(J))
80 CONTINUE
DD=1.0D0
DO 90 J=1,MY
IF(J.NE.I) DD=DD*(DBLE(Y(J))-DAUX)
90 CONTINUE
WORK=((-1.0D0)**(MY-1))*DAUX**(-IY)/DD
DO 101 J=1,MX
DO 100 K=1,MY
PROD2(J,I)=PROD2(J,I)+WORK*DDA(K-1)*PROD1(J,K)
100 CONTINUE
101 CONTINUE
102 CONTINUE
*
DO 125 I=1,MX
DO 120 J=1,MY
OMEG(I,J)=REAL(PROD2(I,J))
120 CONTINUE
125 CONTINUE
RETURN
END
|
