Commit 061207a8 authored by george's avatar george
Browse files

Andrea Pagano's f90 QT library and associated C++ integration and test routines

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2778 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 4d7f60c8
#ifndef QT_H
#define QT_H
#define C_CHAR const char*
#define BLINT int*
#define C_BLINT const int*
#define C_BLDOU const double*
#define BLDOU double*
extern "C"{
void ldm_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT);
void ldsld_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT);
void ldv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void mtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void qt2ld_(BLDOU,C_BLDOU, C_BLINT);
void qt2t_(BLDOU, C_BLDOU, C_BLINT);
void s2d_(BLDOU,C_BLDOU, C_BLINT);
void s2u_(BLDOU,C_BLDOU, C_BLINT);
void td_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tm_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tstt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tu_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tut_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void qtv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void qtsqtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
};
#endif
subroutine LdM(X,Ld,M,n)
! COMPUTATIONAL SUBROUTINE: LdM=Ld*M; Ld lower diagonal matrix; M arbitrary
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: Ld(n,n)
real(8), intent(in) :: M(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*M
do i=2,n
if ((Ld(i,i-1)/=0.0)) then
do j=1,n
X(i,j)=Ld(i,i-1)*M(i-1,j)
end do
end if
end do
return
end subroutine LdM
\ No newline at end of file
subroutine LdSLd(X,Ld,S,n)
! COMPUTATIONAL SUBROUTINE: LdM=Ld*S*Ld; Ld lower diagonal matrix; S symmetric
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: Ld(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4), dimension(1,n) :: vk
integer(4) :: i,j, jj, h
h=0.0
do i=2,n
if ((Ld(i,i-1)/= 0.0)) then
h=h+1.0
vk(1,h)=i
end if
end do
if (h==0.0) then
X=0*S
return
end if
do j=1,h
do jj=1,h
X(vk(1,j),vk(1,jj))=Ld(vk(1,j),vk(1,j)-1)*S(vk(1,j)-1,vk(1,jj)-1)*Ld(vk(1,jj),vk(1,jj)-1)
end do
end do
return
end subroutine LdSLd
subroutine LdV(X,Ld,v,n)
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: ld(n,n)
real(8), intent(in) :: v(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i
X=0*v
do i=2,n
X(i,1)=Ld(i,i-1)*v(i-1,1)
end do
return
end subroutine LdV
subroutine MTt(X,M,Tt,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: M(n,n)
real(8), intent(in) :: Tt(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
X=0.0*M
do i=1,n
do j=1,n
stemp = 0.0
do k = j,n
stemp = stemp + M(i,k+i) * Tt(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine MTt
\ No newline at end of file
subroutine QT2Ld(X,QT,n)
! COMPUTATIONAL SUBROUTINE: extracts lower diagonal from Quasi triangular matrix
! COMPUTATIONAL SUBROUTINE: extracts lower diagonal from Quasi triangular matrix
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i
X=0.0*QT
do i=2,n
if (QT(i,i-1)/=0.0) then
X(i,i-1)=QT(i,i-1)
end if
end do
return
end subroutine QT2Ld
subroutine QT2T(X,QT,n)
! COMPUTATIONAL SUBROUTINE: extracts upper triangular from Quasi triangular matrix
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*QT
do i=1,n
do j=i,n
X(i,j)=QT(i,j)
end do
end do
return
end subroutine QT2T
subroutine QTSQTt(X,QT,S,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*S
do i=1,n
do j=1,n
stemp = 0.0
if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
stemp = QT(i,i-1)*S(i-1,1)*QT(i,i-1)
end if
do h = i,n
do k = j,n
stemp = stemp + QT(i,h) * S(h,k) * QT(j,k)
end do
end do
X(i,j)=stemp
X(j,i)=stemp
end do
end do
return
end subroutine QTSQTt
subroutine QTV(X,QT,V,n)
! COMPUTATIONAL SUBROUTINE: X=QT*V QT upper quasi-triangular; V vector
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(in) :: V(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i,k
real(8) :: stemp
do i=1,n
stemp = 0.0
if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
stemp = QT(i,i-1)*v(i-1,1)
end if
do k = 0,n-i
stemp = stemp + QT(i,k+i) * V(k+i,1)
end do
X(i,1)=stemp
end do
return
end subroutine QTV
subroutine S2D(X,S,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty upper triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i
X=0.0*S
do i=1,n
X(i,i)=sqrt(S(i,i))
end do
return
end subroutine S2D
subroutine S2U(X,S,n)
! COMPUTATIONAL SUBROUTINE: S2U extracts striclty upper triangular from symmetric
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*S
do i=1,n
do j=i+1,n
X(i,j)=S(i,j)
end do
end do
return
end subroutine S2U
subroutine TD(X,T,D,n)
! COMPUTATIONAL SUBROUTINE TUt=T*D T upper triangular; D diagonal
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: D(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*T
do i=1,n
do j=i,n
X(i,j)=T(i,j)*D(j,j)
end do
end do
return
end subroutine TD
subroutine TM(X,T,M,n)
! COMPUTATIONAL SUBROUTINE: TM=T*M T upper triangular; M arbitrary
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: M(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
do i=1,n
do j=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T(i,k+i) * M(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TM
subroutine TSTt(X,T,S,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*S
do i=1,n
do j=1,n
stemp = 0.0
do h = i,n
do k = j,n
stemp = stemp + T(i,h) * S(h,k) * T(j,k)
end do
end do
X(i,j)=stemp
X(j,i)=stemp
end do
end do
return
end subroutine TSTt
subroutine TT(X, T1,T2,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T1(n,n)
real(8), intent(in) :: T2(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*T1
do i=1,n
do j=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T1(i,k+i) * T2(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TT
\ No newline at end of file
subroutine TU(X,T,U,n)
! COMPUTATIONAL SUBROUTINE: TU=T*U; T upper triangular matrix; U strictly upper triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: U(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
X=0.0*T
do i=1,n-1
do j=i+1,n
stemp = 0.0
do k = i,j-1
stemp = stemp + T(i,k) * U(k,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TU
subroutine TUt(X,T,Ut,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: Ut(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*T
do i=1,n
do j=1,n-1
h=max(i,j)
stemp = 0.0
do k = 0,n-h
stemp = stemp + T(i,k+h) * Ut(k+h,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TUt
\ No newline at end of file
subroutine TV(X,T,V,n)
! COMPUTATIONAL SUBROUTINE: TV=T*V T upper triangular; V vector
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: V(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i,k
real(8) :: stemp
do i=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T(i,k+i) * V(k+i,1)
end do
X(i,1)=stemp
end do
return
end subroutine TV
# $Id: Makefile 531 2005-11-30 13:49:48Z kamenik $
# Copyright 2005, Ondra Kamenik
DEBUG = yes
#LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
CC_FLAGS := -DMATLAB -DWINDOWS -DNO_BLAS_H -DNO_LAPACK_H \
-Wall -I../cc -I../sylv/cc -I../cc \
-Ic:/"Program Files"/MATLAB_SV71/extern/include #-pg
ifeq ($(DEBUG),yes)
CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
# CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g
KALMANLIB := kalmanlib_dbg.a
else
# CC_FLAGS := $(CC_FLAGS) -O2
CC_FLAGS := -DTIMING_LOOP $(CC_FLAGS) -O2
KALMANLIB := kalmanlib.a
endif
# Added by GP
# LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
#LDFLAGS := -Wl,--library-path $(LD_LIBRARY_PATH)
LD_LIBS := -Wl,--library-path \
-Wl,-L'f:/MinGW/lib' \
-Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ \
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
-lf95 -lg2c -lmingw32 -lstdc++ $(LDFLAGS) \
-Wl,-L'C:/MinGW/lib/gcc-lib/i686-pc-mingw32/4.0.4' -Wl,-L'C:/MinGW/lib'
# -Wl,-L'f:/CygWin/usr/local/atlas/lib'
# -Wl,-L'f:/CygWin/lib'
# $(LDFLAGS)
# LD_LIBS :=$(LDFLAGS)
# end add
#matrix_interface := GeneralMatrix Vector SylvException
#matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
#mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
#matcppsource := $(patsubst %, ../sylv/cc/%.cpp, $(matrix_interface))
qtf90source := $(wildcard ../f90/*.f90)
qtobjs := $(patsubst %.f90,%.o,$(qtf90source))
cppsource := $(wildcard *.cpp)
hsource := $(wildcard *.h)
objects := $(patsubst %.cpp,%.o,$(cppsource))
dummy.ch:
touch dummy.ch
%.o: %.cpp $(hsource) $(cppsource)
c++ $(CC_FLAGS) -c $*.cpp
qtamvm_exe.exe: qtamvm_exe.o $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o qtamvm_exe.exe qtamvm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
all: $(objects) qtamvm_exe.exe # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
clear:
rm -f *.o
rm -f *.a