TPT' QTSQTt.f90 performance test standalone routines

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2799 ac1d8469-bf42-47a9-8791-bf33cf982152

mex/sources/kalman/qt/f90/QTSQTt.f90

 subroutine QTSQTt(X,QT,S,n)
-!     COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
+!     COMPUTATIONAL SUBROUTINE X=QT*SQT' QT upper quasi-triangular; S symmetric
 
 implicit none
 
@@ -8,31 +8,43 @@ 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
+integer(4) :: i,j,k,h,k0
 
 real(8) :: stemp
 
-
-X=0.0*S
-
-do i=1,n
-do j=1,n
-stemp = 0.0
+X=0*S
+do  i=1,n
+do j=i, n
 if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
-stemp = QT(i,i-1)*S(i-1,1)*QT(i,i-1)
+h=i-1
+else
+h=i
 end if 
-
-do  h = i,n
-do  k = j,n
-stemp = stemp + QT(i,h) * S(h,k) * QT(j,k)
+if (j > 1 .AND. (QT(j,j-1)/= 0.0)) then
+k=j-1
+k0=k 
+else
+k=j
+k0=k
+end if 
+stemp=0
+do while (h <= n)
+do while (k <= n)
+stemp=stemp+QT(i,h)*S(h,k)*QT(j,k)
+k=k+1
 end do
+k=k0
+h=h+1
 end do
 X(i,j)=stemp
+if (i /= j) then
 X(j,i)=stemp
-
+end if 
 end do
-
 end do
 
+
+
+
 return
 end subroutine QTSQTt

mex/sources/kalman/qt/test/Makefile

@@ -6,11 +6,11 @@
 #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 \
+	-Wall -I../cc -I../../sylv/cc -I../cc -mno-cygwin -mwindows \
 	-Ic:/"Program Files"/MATLAB_SV71/extern/include  #-pg
 
 ifeq ($(DEBUG),yes)
-	CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
+	CC_FLAGS := -DDEBUG $(CC_FLAGS) -g -pg
 #	CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g  
 	KALMANLIB := kalmanlib_dbg.a
 else
@@ -25,10 +25,10 @@ endif
 
         LD_LIBS := -Wl,--library-path  \
               -Wl,-L'f:/MinGW/lib' \
+        -Wl,-L'C:/MinGW/lib/gcc-lib/i686-pc-mingw32/4.0.4' -Wl,-L'C:/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' 
+        -lf95 -lg2c -lmingw32  -lstdc++ $(LDFLAGS)
 
 #                  -Wl,-L'f:/CygWin/usr/local/atlas/lib' 
 #                  -Wl,-L'f:/CygWin/lib' 
@@ -53,6 +53,14 @@ dummy.ch:
 %.o: %.cpp $(hsource) $(cppsource)
 	c++ $(CC_FLAGS) -c $*.cpp
 
+dgemmmtm_exe.exe: dgemmmtm_exe.o $(objects) $(qtobjs) $(hsource) $(cppsource) 
+	gcc $(CC_FLAGS) -o dgemmmtm_exe.exe dgemmmtm_exe.o ascii_array.o \
+	$(qtobjs) $(LD_LIBS) 
+
+qtmmmtm_exe.exe: qtmmmtm_exe.o $(objects) $(qtobjs) $(hsource) $(cppsource) 
+	gcc $(CC_FLAGS) -o qtmmmtm_exe.exe qtmmmtm_exe.o ascii_array.o \
+	$(qtobjs) $(LD_LIBS) 
+
 dgemvm_exe.exe: dgemvm_exe.o $(qtobjs) $(hsource) $(cppsource) 
 	gcc $(CC_FLAGS) -o dgemvm_exe.exe dgemvm_exe.o ascii_array.o \
 	$(qtobjs) $(LD_LIBS) 
@@ -69,7 +77,7 @@ qtv1mvm_exe.exe: qtv1mvm_exe.o $(qtobjs) $(hsource) $(cppsource)
 	gcc $(CC_FLAGS) -o qtv1mvm_exe.exe qtv1mvm_exe.o ascii_array.o \
 	$(qtobjs) $(LD_LIBS) 
 
-all: $(objects) qtv1mvm_exe.exe dgemvm_exe.exe qtvmvm_exe.exe qtamvm_exe.exe # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
+all: $(objects) dgemmmtm_exe.exe qtmmmtm_exe.exe qtv1mvm_exe.exe dgemvm_exe.exe qtvmvm_exe.exe qtamvm_exe.exe # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
 
 clear:
 	rm -f *.o

mex/sources/kalman/qt/test/dgemmmtm_exe.cpp

+/*
+* Copyright (C) 2008-2009 Dynare Team
+*
+* This file is part of Dynare.
+*
+* Dynare is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* Dynare is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/****************************************************** 
+%
+% This provides an interface to BLAS dgemv f90 library function  
+% to multiply Quasi trinagular matrix (T) with a vector a
+%
+%
+% use:
+%     dgemmmtm_exe QTt_file SS_file size [loops - if enabled]
+%
+% NOTE: due to fortran matrix orientation, input matrices need to be passed 
+% as transposed so QTt instead QT
+%
+%
+% INPUTS
+%    QT                      [double]    mm*mm transition matrix of the state equation.
+%    SS                      [double]    mm*mm state cov matrix.
+%
+% OUTPUTS
+%    TSTt update             [double]    mm*mm state cov matrix updated.
+%                                       as file: a_file_out
+**********************************************************/
+
+#include "cppblas.h"
+#include <stdio.h>
+#include <math.h>
+#include "ascii_array.h"
+#include <iostream>
+#include <stdexcept>
+#include <malloc.h>
+
+int main(int argc, char* argv[])
+    {
+    
+    if (argc < 3 )
+      {
+      printf("Must have  min 2 input parameters.\n");
+      exit(1);
+      }
+  
+    try 
+      {
+      // make input matrices
+      int n=atoi(argv[3]);
+      double *TSTt, *TS ;
+      AsciiNumberArray QT, SS;
+      QT.GetMX(argv[1],n,n);
+      SS.GetMX(argv[2],n,n);
+      const double alpha=1.0;
+      const double beta=0.0;
+      // create output and upload output data
+      TS=(double *)calloc(n*n, sizeof(double));
+      TSTt=(double *)calloc(n*n, sizeof(double));
+
+      
+#ifdef TIMING_LOOP
+      int loops=1;//000;
+      if (argc >3 )
+        loops = atoi(argv[4]);
+      for (int tt=0;tt<loops;++tt)
+        {
+#endif
+
+/* DGEMM  performs one of the matrix-matrix operations
+*
+*     C := alpha*op( A )*op( B ) + beta*C,
+*
+*  where  op( X ) is one of
+*
+*     op( X ) = X   or   op( X ) = X',
+*
+* void BLAS_dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n,
+*					CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
+*					CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
+*					BLDOU c, CONST_BLINT ldc);
+*
+* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+*
+*/
+	  BLAS_dgemm("N", "N",  &n,  &n, &n, &alpha, QT.data, &n, SS.data, &n, &beta, TS, &n);
+
+	  BLAS_dgemm("N", "T", &n,  &n, &n, &alpha, TS, &n, QT.data, &n, &beta, TSTt, &n);
+
+#ifdef TIMING_LOOP
+        }
+      printf("QT array mvm: finished: %d loops\n",loops);
+#endif
+      // create output and upload output data
+      WriteMX(argv[2], TSTt,n,n);
+      free(TSTt);
+      free(TS);
+      } 
+    catch (std::exception e) 
+      {
+      std::cout <<"Error" << std::endl;
+      }
+    
+  }; //main

mex/sources/kalman/qt/test/qtmmmtm_exe.cpp

+/*
+* Copyright (C) 2008-2009 Dynare Team
+*
+* This file is part of Dynare.
+*
+* Dynare is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* Dynare is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/****************************************************** 
+%
+% This provides an interface to BLAS dgemm  library function  
+% to multiply Quasi trinagular matrix (T) with another matrix
+%
+% use:
+%     dgemmmtm_exe QTt_file SS_file size [loops - if enabled]
+%
+% NOTE: due to fortran matrix orientation, input matrices need to be passed 
+% as transposed so QTt instead QT
+%
+%
+% INPUTS
+%    QT                      [double]    mm*mm transition matrix of the state equation.
+%    SS                      [double]    mm*mm state cov matrix.
+%
+% OUTPUTS
+%    TSTt update             [double]    mm*mm state cov matrix updated.
+%                                       as file: a_file_out
+**********************************************************/
+
+#include "qt.h"
+#include <stdio.h>
+#include <math.h>
+#include "ascii_array.h"
+#include <iostream>
+#include <stdexcept>
+#include <malloc.h>
+
+int main(int argc, char* argv[])
+    {
+    
+    if (argc < 3 )
+      {
+      printf("Must have  min 2 input parameters.\n");
+      exit(1);
+      }
+  
+    try 
+      {
+      // make input matrices
+      int n=atoi(argv[3]);
+      double *TSTt ;
+      AsciiNumberArray QT, SS;
+      QT.GetMX(argv[1],n,n);
+      SS.GetMX(argv[2],n,n);
+      // create output and upload output data
+      TSTt=(double *)calloc(n*n, sizeof(double));
+
+      
+#ifdef TIMING_LOOP
+      int loops=1;//000;
+      if (argc >3 )
+        loops = atoi(argv[4]);
+      for (int tt=0;tt<loops;++tt)
+        {
+#endif
+
+/* qtsqtt_  performs one of the matrix-matrix operations
+*
+*     C := QT*S*QT'
+*/
+
+      qtsqtt_(TSTt, QT.data, SS.data, &n);
+
+
+#ifdef TIMING_LOOP
+        }
+      printf("QT array mvm: finished: %d loops\n",loops);
+#endif
+      // create output and upload output data
+      WriteMX(argv[2], TSTt,n,n);
+      free(TSTt);
+      } 
+    catch (std::exception e) 
+      {
+      std::cout <<"Error" << std::endl;
+      }
+    
+  }; //main