Refactoring kalman.cpp to use dgemm more directly and bypass some of...

Refactoring kalman.cpp to use dgemm more directly and bypass some of re-casting matrix constructors.

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

mex/sources/kalman/cc/Makefile

-# $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../../gensylv/cc  \
-		-Ic:/"Program Files"/MATLAB_SV71/extern/include
-
-ifeq ($(DEBUG),yes)
-	CC_FLAGS :=  -DDEBUG $(CC_FLAGS) -g
-else
-	CC_FLAGS := $(CC_FLAGS) -O2
-endif
-
-# Added by GP
-#    LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
-        LDFLAGS := -Wl,--library-path $(LD_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 \
-        -lg2c -lmingw32  -lstdc++ $(LDFLAGS)
-
-#                  -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 %, ../../gensylv/cc/%.o, $(matrix_interface))
-mathsource := $(patsubst %,  ../../gensylv/cc/%.h, $(matrix_interface))
-matcppsource := $(patsubst %,  ../../gensylv/cc/%.cpp, $(matrix_interface))
-cwebsource := $(wildcard *.cweb)
-# cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource)) 
-cppsource := $(wildcard *.cpp)
-hsource := $(wildcard *.h)
-objects := $(patsubst %.cpp,%.o,$(cppsource))
-hwebsource := $(wildcard *.hweb)
-
-dummy.ch:
-	touch dummy.ch
-
-# %.cpp: %.cweb dummy.ch
-# 	ctangle -bhp $*.cweb dummy.ch $*.cpp
-
-# %.h: %.hweb dummy.ch
-# 	ctangle -bhp $*.hweb dummy.ch $*.h
-
-%.o: %.cpp $(hsource) $(mathsource)
-	c++ $(CC_FLAGS) -c $*.cpp
-
-all: $(objects) # $(cppsource) $(hsource)
-
-doc: main.web $(hwebsource) $(cwebsource)
-	cweave -bhp main.web
-	pdftex main
-	mv main.pdf ts.pdf 
-
-clear:
-	rm -f *.o
-	rm -f *.{pdf,dvi,log,scn,idx,toc}
-#	rm -f *.cpp
-#	rm -f *.h
+# $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../sylv/cc  \
+		-Ic:/"Program Files"/MATLAB_SV71/extern/include #-pg
+
+ifeq ($(DEBUG),yes)
+	CC_FLAGS :=  -DDEBUG $(CC_FLAGS) -g
+else
+	CC_FLAGS := $(CC_FLAGS) -O3
+endif
+
+# Added by GP
+#    LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
+        LDFLAGS := -Wl,--library-path $(LD_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 \
+        -lg2c -lmingw32  -lstdc++ $(LDFLAGS)
+
+#                  -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))
+cwebsource := $(wildcard *.cweb)
+# cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource)) 
+cppsource := $(wildcard *.cpp)
+hsource := $(wildcard *.h)
+objects := $(patsubst %.cpp,%.o,$(cppsource))
+hwebsource := $(wildcard *.hweb)
+
+dummy.ch:
+	touch dummy.ch
+
+# %.cpp: %.cweb dummy.ch
+# 	ctangle -bhp $*.cweb dummy.ch $*.cpp
+
+# %.h: %.hweb dummy.ch
+# 	ctangle -bhp $*.hweb dummy.ch $*.h
+
+%.o: %.cpp $(hsource) $(mathsource)
+	c++ $(CC_FLAGS) -c $*.cpp
+
+all: $(objects) # $(cppsource) $(hsource)
+
+doc: main.web $(hwebsource) $(cwebsource)
+	cweave -bhp main.web
+	pdftex main
+	mv main.pdf ts.pdf 
+
+clear:
+	rm -f *.o
+	rm -f *.{pdf,dvi,log,scn,idx,toc}
+#	rm -f *.cpp
+#	rm -f *.h

mex/sources/kalman/cc/kalman.cpp

@@ -29,6 +29,9 @@ routines.
 #include "kalman.h"
 #include "ts_exception.h"
 
+#include "cppblas.h"
+//#include "cpplapack.h"
+
 #include <math.h> 
 #include <float.h> 
 #include <cmath> 
@@ -609,7 +612,8 @@ BasicKalmanTask::filter(int&per,int&d, int start, std::vector<double>* vll)const
   per= vll->size() ;
   return filterNonDiffuse(init.getA(),init.getPstar(), start, vll);
   }
-  /*****************
+
+/*****************
   This is public interface that runs a filter followed by a smoother. 
   In addition to things returned by |KalmanTask::filter|, it fills also 
   |SmootherResults|, which must be initialized to the number of columns 
@@ -648,7 +652,7 @@ KalmanTask::filter_and_smooth(SmootherResults&sres,
   }
 
 
-  /*****************
+/*****************
   This runs a Basic non-diffuse filter with the given $t$, $a_t$ and
   $P_t$. It fills the |FilterResults|.
   
@@ -670,24 +674,24 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
   double loglik=0;
   Vector at(a);
   GeneralMatrix Pt(P);
-  GeneralMatrix PtZeros(Pt.numRows(), Pt.numCols());
-  PtZeros.zeros();
+//  GeneralMatrix PtZeros(Pt.numRows(), Pt.numCols());
+//  PtZeros.zeros();
   if(TSUtils::hasNegativeDiagonal(Pt)||!TSUtils::isSymDiagDominant(Pt))
     return 0.0;
-    /*
-    ConstGeneralMatrix Zt(Z);
-    ConstGeneralMatrix Ht(H);
-    ConstGeneralMatrix Tt(T);
-    ConstGeneralMatrix Qt(Q);
-    ConstGeneralMatrix Rt(R);
-  */
   GeneralMatrix Ft (Ht.numRows(), Ht.numCols()  );
+  GeneralMatrix Lt(Tt);
 //  PLUFact Ftinv(Ft);
 
   bool isTunit=0;// Tt->isUnit();
   bool isQzero= Qt.isZero();
   bool isRzero= Rt.isZero();
-  
+  const double alpha=1.0;
+  const double neg_alpha=-1.0;
+  const double omega=0.0;
+  const int m=Pt.numRows();
+  const int n=Zt.numRows();
+  const int rcols= Rt.numCols();
+
   int t= 1; 
   int nonsteady=1;
   for(;t<=data.numCols()&&nonsteady;++t)
@@ -706,7 +710,10 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
     GeneralMatrix Mt(Pt,Zt,"trans");
 //    GeneralMatrix Ft(Ht);
     Ft=Ht;
-    Ft.multAndAdd(Zt,ConstGeneralMatrix(Mt));
+//    Ft.multAndAdd(Zt,ConstGeneralMatrix(Mt));
+    // DGEMM: C := alpha*op( A )*op( B ) + beta*C,
+		BLAS_dgemm("N", "N", &n, &n, &m, &alpha, Zt.base(), &n,
+				   Mt.base(), &m, &alpha, Ft.base(), &n); 
     
     PLUFact Ftinv(Ft); //    Ftinv=Ft;
     if(!Ftinv.isRegular())
@@ -721,8 +728,12 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
     /*****************
     This calculates  $$L_t = T_t-K_tZ_t.$$
     *****************/   
-    GeneralMatrix Lt(Tt);
-    Lt.multAndAdd(ConstGeneralMatrix(Kt),Zt,-1.0);
+    //GeneralMatrix Lt(Tt);
+    Lt=Tt;
+    //Lt.multAndAdd(ConstGeneralMatrix(Kt),Zt,-1.0);
+    // DGEMM: C := alpha*op( A )*op( B ) + beta*C,
+		BLAS_dgemm("N", "N", &m, &m, &n, &neg_alpha, Kt.base(), &m,
+				   Zt.base(), &n, &alpha, Lt.base(), &m); 
     
     
     /*****************
@@ -752,9 +763,11 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
       if(!isTunit)
         {
 //        Pt.zeros();
-        Pt=PtZeros;
-        Pt.multAndAdd(Tt,ConstGeneralMatrix(PtLttrans));
-//        Pt=mult(Tt,ConstGeneralMatrix(PtLttrans));        
+//        Pt=PtZeros;
+//        Pt.multAndAdd(Tt,ConstGeneralMatrix(PtLttrans));
+        // DGEMM: C := alpha*op( A )*op( B ) + beta*C,
+		    BLAS_dgemm("N", "N", &m, &m, &m, &alpha, Tt.base(), &m,
+				       PtLttrans.base(), &m, &omega, Pt.base(), &m); 
         }
       else
         {
@@ -763,7 +776,10 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
       if(!isRzero&&!isQzero)
         {
         GeneralMatrix QtRttrans(Qt,Rt,"trans");
-        Pt.multAndAdd(Rt,ConstGeneralMatrix(QtRttrans));
+ //       Pt.multAndAdd(Rt,ConstGeneralMatrix(QtRttrans));
+        // DGEMM: C := alpha*op( A )*op( B ) + beta*C,
+		    BLAS_dgemm("N", "N", &m, &m,  &rcols, &alpha, Rt.base(), &m,
+				       QtRttrans.base(), &rcols, &alpha, Pt.base(), &m); 
         }
       
       }

mex/sources/kalman/matlab/Makefile

@@ -6,32 +6,34 @@
 #LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
 
 CC_FLAGS :=  -DMATLAB -DWINDOWS  -DNO_BLAS_H -DNO_LAPACK_H \
-	-Wall -I../qt/cc -I../sylv/cc -I../cc \
+	-Wall -mno-cygwin -I../qt/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 #-pg #-Wl,-pg
+	CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
+#	CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g #-pg #-Wl,-pg
 	KALMANLIB := kalmanlib_dbg.a
 else
-#	CC_FLAGS := $(CC_FLAGS) -O2
-	CC_FLAGS := -DTIMING_LOOP $(CC_FLAGS) -O2
+#	CC_FLAGS := $(CC_FLAGS) -O3
+	CC_FLAGS := -DTIMING_LOOP $(CC_FLAGS) -O3
 	KALMANLIB := kalmanlib.a
 endif
 
 # Added by GP
 #    LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
-        #LDFLAGS := -Wl,--library-path  $(LD_LIBRARY_PATH) 
+        #LDFLAG := -Wl,--library-path  $(LD_LIBRARY_PATH) 
 
         LD_LIBS := -Wl,--library-path  \
-              -Wl,-L'f:/MinGW/lib' \
+                  -Wl,-L'f:/CygWin/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'C:/MinGW/lib/gcc-lib/i686-pc-mingw32/4.0.4' \
+	-Wl,-L'C:/MinGW/lib'  kalmanlib.def
 
 #                  -Wl,-L'f:/CygWin/usr/local/atlas/lib' 
 #                  -Wl,-L'f:/CygWin/lib' 
+#              -Wl,-L'f:/MinGW/lib' 
 #                 $(LDFLAGS)
 #    LD_LIBS :=$(LDFLAGS)
 # end add
@@ -93,7 +95,7 @@ qtmvm.dll: qtmvm.o $(KALMANLIB) # $(hsource) $(cppsource)
 
 kalman_filter_dll.dll: kalman_filters.o $(KALMANLIB) # $(hsource) $(cppsource) 
 	gcc -shared $(CC_FLAGS) -o kalman_filter_dll.dll kalman_filters.o \
-	$(KALMANLIB) $(LD_LIBS) 
+	$(KALMANLIB) $(LD_LIBS)
 
 kalman_filters_testx.exe: kalman_filters_testx.o $(KALMANLIB) # $(hsource) $(cppsource) 
 	gcc $(CC_FLAGS) -pg -o kalman_filters_testx.exe kalman_filters_testx.o \

mex/sources/kalman/matlab/kalman_filters_testx.cpp

@@ -58,6 +58,9 @@
 % NOTES
 %   The vector "vll" is used to evaluate the jacobian of the likelihood.
 **********************************************************/
+
+#include <iostream>
+using namespace std;
 
 #include "kalman.h"
 #include "ts_exception.h"
@@ -65,6 +68,7 @@
 #include "GeneralMatrix.h"
 #include "Vector.h"
 #include "SylvException.h"
+
 
 #include "mex.h"
 
@@ -263,10 +267,11 @@ int main(int argc, char* argv[])
 #ifdef TIMING_LOOP
   for (int tt=0;tt<10000;++tt)
     {
-#endif
+#endif 
         loglik = bkt.filter( per, d, (start-1), vll);
 #ifdef DEBUG		
-    mexPrintf("Basickalman_filter: loglik=%d \n", loglik);
+//    mexPrintf("Basickalman_filter: loglik=%f \n", loglik);
+//    cout << "loglik " << loglik << "\n";
 #endif		
 #ifdef TIMING_LOOP
     }