Initial commit from Polytechnique Montreal

1 files changed, 119 insertions, 0 deletions

diff --git a/Dragon/src/LIBOMG.f b/Dragon/src/LIBOMG.f
new file mode 100644
index 0000000..e4b7162
--- /dev/null
+++ b/Dragon/src/LIBOMG.f
@@ -0,0 +1,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