Added multInvRight function and temporary matrix space for reducing...

Added multInvRight function and temporary matrix space for reducing constructor calls to GeneralMatrix

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

mex/sources/kalman/sylv/cc/GeneralMatrix.cpp

@@ -85,6 +85,8 @@ GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum1,
 
 GeneralMatrix::~GeneralMatrix()
   {
+     if(tmpGMp)
+       delete tmpGMp;
   }
 
 
@@ -172,6 +174,52 @@ void GeneralMatrix::multLeftTrans(const ConstGeneralMatrix& m)
   gemm_partial_left("T", m, 1.0, 0.0);
   }
 
+
+/* this = this * A^(-1)  */
+void 
+GeneralMatrix::multInvRight( GeneralMatrix&A)
+  {
+  // check or allocate tmp space for Transpose *this
+  if (tmpGMp)
+    {
+    if (tmpGMp->numCols()!=rows || tmpGMp->numRows()!=cols)
+      delete (tmpGMp);
+    }
+  if (!tmpGMp) 
+    tmpGMp= new GeneralMatrix(cols,rows); // allocate space only once if and when needed!
+
+  // tmpGMp=(*this)' i.e. Transpose (*this)
+  for (int i = 0; i < rows; i++)
+    for (int j = 0; j < cols; j++)
+      tmpGMp->get(j,i) = get(i,j);
+  // check A and this suiability
+  const int mcols=A.numCols();
+  if (A.numRows() != mcols) {
+    throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
+    }
+  if (cols != A.numRows()) {
+    throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
+    }
+  
+  if (rows > 0) 
+    {
+    /* out =tmpGMp out = inv(A')*(*this)' = inv(A')*(*tmpGMp) */
+    int* ipiv = new int[cols];
+    int info;
+    const int mld=A.getLD();
+    LAPACK_dgetrf(&cols, &cols, A.base(), &cols, ipiv, &info);
+    LAPACK_dgetrs("T", &cols, &mcols, A.base(), &cols, ipiv, tmpGMp->base(),
+      &mld, &info);
+    delete [] ipiv;
+    // *this= Transpose(tmpGMp out)
+    for (int i = 0; i < rows; i++)
+      for (int j = 0; j < cols; j++)
+        get(i,j) = tmpGMp->get(j,i);
+    }
+
+  }
+
+
 // here we must be careful for ld
 void GeneralMatrix::zeros()
   {

mex/sources/kalman/sylv/cc/GeneralMatrix.h

@@ -76,7 +76,7 @@ class GeneralMatrix {
     int rows;
     int cols;
     int ld;
-
+    GeneralMatrix * tmpGMp;
 public:
     GeneralMatrix(int m, int n)
       : data(m*n), rows(m), cols(n), ld(m) {}
@@ -195,6 +195,9 @@ public:
     void multLeftTrans(const GeneralMatrix& m)
       {multLeftTrans(ConstGeneralMatrix(m));}
 
+    /* this = this * m^(-1)  */
+    void multInvRight(GeneralMatrix&m);
+    
     /* x = scalar(a)*x + scalar(b)*this*d */
     void multVec(double a, Vector& x, double b, const ConstVector& d) const;
 //      {ConstGeneralMatrix(*this).multVec(a, x, b, d);}

mex/sources/kalman/sylv/cc/Makefile

@@ -2,8 +2,8 @@
 
 # Tag $Name:  $
 
-CC_FLAGS := -Wall -I../testing -I../cc -DMATLAB  \
-        $(CC_INCLUDE_PATH) -Ic:/"Program Files"/MATLAB_SV71/extern/include
+CC_FLAGS := -Wall -Winline -I../testing -I../cc -DMATLAB  \
+        $(CC_INCLUDE_PATH) -Ic:/"Program Files"/MATLAB_SV71/extern/include #-pg
 
 LDFLAGS = -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ \
 	-Wl,-llibmex -Wl,-llibmx  -Wl,-llibmwlapack -Wl,-llibdflapack \
@@ -15,7 +15,7 @@ ifeq ($(DEBUG),yes)
 #	CC_FLAGS := $(CC_FLAGS) -g -DTL_DEBUG=2
     CC_FLAGS := $(CC_FLAGS) -g -DPOSIX_THREADS
 else
-	CC_FLAGS := $(CC_FLAGS) -O2
+	CC_FLAGS := $(CC_FLAGS) -O3
 endif
 
 ifeq ($(OS),Windows_NT)