Updates to Kalman filter library and dll and dll tests

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

mex/sources/kalman/cc/kalman.cpp

@@ -670,6 +670,8 @@ 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();
   if(TSUtils::hasNegativeDiagonal(Pt)||!TSUtils::isSymDiagDominant(Pt))
     return 0.0;
     /*
@@ -679,6 +681,9 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
     ConstGeneralMatrix Qt(Q);
     ConstGeneralMatrix Rt(R);
   */
+  GeneralMatrix Ft (Ht.numRows(), Ht.numCols()  );
+//  PLUFact Ftinv(Ft);
+
   bool isTunit=0;// Tt->isUnit();
   bool isQzero= Qt.isZero();
   bool isRzero= Rt.isZero();
@@ -699,11 +704,11 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
     This calculates $$F_t = Z_tP_tZ_t^T+H_t.$$
     *****************/   
     GeneralMatrix Mt(Pt,Zt,"trans");
-    GeneralMatrix Ft(Ht);
+//    GeneralMatrix Ft(Ht);
+    Ft=Ht;
     Ft.multAndAdd(Zt,ConstGeneralMatrix(Mt));
     
-    
-    PLUFact Ftinv(Ft);
+    PLUFact Ftinv(Ft); //    Ftinv=Ft;
     if(!Ftinv.isRegular())
       return 0.0;
     
@@ -746,8 +751,10 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
       GeneralMatrix PtLttrans(Pt,Lt,"trans");
       if(!isTunit)
         {
-        Pt.zeros();
+//        Pt.zeros();
+        Pt=PtZeros;
         Pt.multAndAdd(Tt,ConstGeneralMatrix(PtLttrans));
+//        Pt=mult(Tt,ConstGeneralMatrix(PtLttrans));        
         }
       else
         {

mex/sources/kalman/matlab/Makefile

@@ -4,9 +4,10 @@
 #DEBUG = yes
 
 #LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
-CC_FLAGS :=  -DMATLAB -DWINDOWS  -DNO_BLAS_H -DNO_LAPACK_H  \
-	-Wall -I../sylv/cc -I../cc \
-	-Ic:/"Program Files"/MATLAB_SV71/extern/include 
+
+CC_FLAGS :=  -DMATLAB -DWINDOWS  -DNO_BLAS_H -DNO_LAPACK_H \
+	-Wall -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
@@ -20,15 +21,17 @@ endif
 
 # Added by GP
 #    LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
-        #LDFLAGS := -Wl,--library-path 
-        LD_LIBS := -Wl,--library-path $(LD_LIBRARY_PATH) \
-             -Wl,-L'f:/MinGW/lib' \
+        #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 \
-        -lg2c -lmingw32  -lstdc++ $(LDFLAGS)
+        -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' \
+#                  -Wl,-L'f:/CygWin/usr/local/atlas/lib' 
+#                  -Wl,-L'f:/CygWin/lib' 
 #                 $(LDFLAGS)
 #    LD_LIBS :=$(LDFLAGS)
 # end add
@@ -37,6 +40,8 @@ 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 ../qt/f90/*.f90)
+qtobjs := $(patsubst %.f90,%.o,$(qtf90source))
 # cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource)) 
 kalmancppsource := $(wildcard ../cc/*.cpp)
 kalmanhsource := $(wildcard ../cc/*.h)
@@ -62,8 +67,8 @@ dummy.ch:
 	c++ $(CC_FLAGS) -c $*.cpp
 
 
-$(KALMANLIB): $(objects) # $(matobjs) #$(kalmanobjects)
-	ar cr $(KALMANLIB) $(kalmanobjects) $(matobjs) # $(objects)  
+$(KALMANLIB): $(objects) # $(matobjs) $(qtobjs) #$(kalmanobjects)
+	ar cr $(KALMANLIB) $(kalmanobjects) $(matobjs) $(qtobjs)  
 	ranlib $(KALMANLIB)
 
 kalman_smoother_dll.dll: kalman_smoother.o $(KALMANLIB) #$(hsource) $(cppsource) 
@@ -74,6 +79,18 @@ minv.dll: minv.o $(KALMANLIB) # $(hsource) $(cppsource)
 	gcc -shared $(CC_FLAGS) -o minv.dll minv.o \
 	kalmanlib.a $(LD_LIBS) 
 
+gmvm.dll: gmvm.o $(KALMANLIB) # $(hsource) $(cppsource) 
+	gcc -shared $(CC_FLAGS) -o gmvm.dll gmvm.o \
+	kalmanlib.a $(LD_LIBS) 
+
+qtamvm.dll: qtamvm.o $(KALMANLIB) # $(hsource) $(cppsource) 
+	gcc -shared $(CC_FLAGS) -o qtamvm.dll qtamvm.o \
+	kalmanlib.a $(LD_LIBS) 
+
+qtmvm.dll: qtmvm.o $(KALMANLIB) # $(hsource) $(cppsource) 
+	gcc -shared $(CC_FLAGS) -o qtmvm.dll qtmvm.o \
+	kalmanlib.a $(LD_LIBS) 
+
 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) 

mex/sources/kalman/matlab/gmvm.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 QT f90 library by Andrea Pagano 
+% to multiply Quasi trinagular matrix (T) with a vector a
+%
+% function [a] = qtmvm(QT,a)
+%
+%  1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
+%  2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+%
+% OUTPUTS
+%    Tinverse               [double]    mm*mm transition matrix of the state equation.
+**********************************************************/
+
+#include "qt.h"
+#include "kalman.h"
+#include "ts_exception.h"
+
+#include "GeneralMatrix.h"
+#include "Vector.h"
+#include "SylvException.h"
+
+#include "mex.h"
+
+extern "C" {
+  void mexFunction(int nlhs, mxArray* plhs[],
+    int nrhs, const mxArray* prhs[])
+    {
+    
+    if (nrhs < 2 )
+      mexErrMsgTxt("Must have  min 2 input parameters.\n");
+    if (nlhs < 1 )
+      mexErrMsgTxt("Must have  min 1 output parameters.\n")
+      
+      ;
+    //int start = 1; // default start of likelihood calculation
+    // test for univariate case
+    try 
+      {
+      // make input matrices
+      int n=mxGetM(prhs[0]);
+      
+      ConstGeneralMatrix QT(mxGetPr(prhs[0]), n, mxGetN(prhs[0]));
+      //      ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      ConstVector a (mxGetPr(prhs[1]), n);
+      
+      // create output and upload output data
+      plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
+      //      double * mxinv= mxGetPr(plhs[0]);
+      //      GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      Vector Ta (mxGetPr(plhs[0]), n);
+      //    Tinv.unit();
+      //      Ta.zeros();
+      
+      
+      // make storage for output
+      
+#ifdef TIMING_LOOP
+      int loops=1;//000;
+      if (nrhs >2 )
+        loops = (int)mxGetScalar(prhs[2]);
+      for (int tt=0;tt<loops;++tt)
+        {
+#endif
+#ifdef DEBUG	
+      //  QT.print();
+#endif
+        QT.multVec(0.0, Ta, 1.0, a);
+
+#ifdef DEBUG	
+        Ta.print();
+#endif        
+#ifdef TIMING_LOOP
+        }
+      mexPrintf("gmvm: finished: %d loops\n",loops);
+#endif
+      // create output and upload output data
+      /*      if (nlhs >= 1)
+      {
+      plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
+      double * mxinv= mxGetPr(plhs[0]);
+      // allocate likelihood array
+      for (int j=0;j<nper;++j)
+      mxinv[j]=(*vll)[j];
+      }
+      */      
+      } 
+    catch (const TSException& e) 
+      {
+      mexErrMsgTxt(e.getMessage());
+      } 
+    catch (SylvException& e) 
+      {
+      char mes[300];
+      e.printMessage(mes, 299);
+      mexErrMsgTxt(mes);
+      }
+    
+    } // mexFunction
+  }; // extern 'C'

mex/sources/kalman/matlab/kalman_filters.cpp

@@ -190,6 +190,7 @@ extern "C" {
         for (int j=0;j<nper;++j)
           mxll[j]=(*vll)[j];
         }
+      delete vll;
       
       } 
     catch (const TSException& e) 

mex/sources/kalman/matlab/kalman_filters_testx.cpp

@@ -278,7 +278,8 @@ int main(int argc, char* argv[])
 
         // destroy init
       delete init;
-      mexPrintf("logLik = %f \n", loglik);
+      mexPrintf("logLik = %f \n", loglik);
+      delete vll;
       // create output and upload output data
 /************
       if (nlhs >= 1)

mex/sources/kalman/matlab/minv.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/>.
+*/
+
+/* derived from c++kalman_filter library by O. Kamenik */
+
+// This provides a matrix inversion
+
+/****************************************************** 
+% function [Tinv] = minv(T)
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+
+% OUTPUTS
+%    Tinverse               [double]    mm*mm transition matrix of the state equation.
+**********************************************************/
+
+#include "kalman.h"
+#include "ts_exception.h"
+
+#include "GeneralMatrix.h"
+#include "Vector.h"
+#include "SylvException.h"
+
+#include "mex.h"
+
+extern "C" {
+  void mexFunction(int nlhs, mxArray* plhs[],
+    int nrhs, const mxArray* prhs[])
+    {
+
+    if (nrhs < 1 )
+      mexErrMsgTxt("Must have  min 1 input parameters.\n");
+    if (nlhs < 1 )
+      mexErrMsgTxt("Must have  min 1 output parameters.\n")
+      
+      ;
+    //int start = 1; // default start of likelihood calculation
+    // test for univariate case
+    try 
+      {
+      // make input matrices
+      ConstGeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
+      // create output and upload output data
+      plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxREAL);
+//      double * mxinv= mxGetPr(plhs[0]);
+      GeneralMatrix Tinv(mxGetPr(plhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
+//    Tinv.unit();
+//    Tinv.zeros();
+
+
+      // make storage for output
+
+#ifdef TIMING_LOOP
+      int loops=1000;
+      if (nrhs >1 )
+         loops = (int)mxGetScalar(prhs[1]);
+      for (int tt=0;tt<loops;++tt)
+      {
+#endif
+        Tinv.unit();
+        T.multInvLeft(Tinv);
+        //Tinv.print();
+
+#ifdef TIMING_LOOP
+      }
+      mexPrintf("minv: finished: %d loops\n",loops);
+#endif
+      // create output and upload output data
+/*      if (nlhs >= 1)
+        {
+        plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
+        double * mxinv= mxGetPr(plhs[0]);
+        // allocate likelihood array
+        for (int j=0;j<nper;++j)
+          mxinv[j]=(*vll)[j];
+        }
+*/      
+      } 
+    catch (const TSException& e) 
+      {
+      mexErrMsgTxt(e.getMessage());
+      } 
+    catch (SylvException& e) 
+      {
+      char mes[300];
+      e.printMessage(mes, 299);
+      mexErrMsgTxt(mes);
+      }
+
+    } // mexFunction
+  }; // extern 'C'

mex/sources/kalman/matlab/qtamvm.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 QT f90 library by Andrea Pagano 
+% to multiply Quasi trinagular matrix (T) with a vector a
+%
+% function [a] = qtmvm(QT,a)
+%
+%  1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
+%  2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+%    a                      [double]    mm state vector.
+%
+% OUTPUTS
+%    Tinverse               [double]    mm*mm transition matrix of the state equation.
+**********************************************************/
+
+#include "qt.h"
+#include "kalman.h"
+#include "ts_exception.h"
+
+#include "GeneralMatrix.h"
+#include "Vector.h"
+#include "SylvException.h"
+
+#include "mex.h"
+
+extern "C" {
+  void mexFunction(int nlhs, mxArray* plhs[],
+    int nrhs, const mxArray* prhs[])
+    {
+    
+    if (nrhs < 2 )
+      mexErrMsgTxt("Must have  min 2 input parameters.\n");
+    if (nlhs < 1 )
+      mexErrMsgTxt("Must have  min 1 output parameters.\n")
+      
+      ;
+    //int start = 1; // default start of likelihood calculation
+    // test for univariate case
+    try 
+      {
+      // make input matrices
+      int n=mxGetM(prhs[0]);
+
+      double *T1, *Ld, *TV ;     
+      double * QT=mxGetPr(prhs[0]);
+      //      ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      double *a =(mxGetPr(prhs[1]));
+      
+      // create output and upload output data
+      plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
+      //      double * mxinv= mxGetPr(plhs[0]);
+      //      GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      double * Ta =mxGetPr(plhs[0]);
+      //    Tinv.unit();
+      //      Ta.zeros();
+
+      T1=(double *)mxCalloc(n*n, sizeof(double));
+      Ld=(double *)mxCalloc(n*n,sizeof(double));
+      TV=(double *)mxCalloc(n, sizeof(double));
+
+      
+#ifdef TIMING_LOOP
+      int loops=1;//000;
+      if (nrhs >2 )
+        loops = (int)mxGetScalar(prhs[2]);
+      for (int tt=0;tt<loops;++tt)
+        {
+#endif
+#ifdef DEBUG	
+      //  QT.print();
+#endif
+        // 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
+     //   double *T1, *Ld, *dTa;//, dT1=-7.77;
+//  mexPrintf("start dT1 = %f\n", dT1);
+//        dT1 = qt2t_(QT.base() ,&n) ;
+// mexPrintf("end dT1 = %f\n", dT1);
+        //T1=&dT1;
+
+        qt2t_(T1, QT ,&n) ;
+//        T1=*T1p;
+//        GeneralMatrix T1gm(T1,n,n);
+//        Ld = qt2ld_(QT.base(),&n);
+        qt2ld_(Ld , QT,&n);
+        // 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
+//        dTa = ldv_(Ld,a.base() ,&n);
+        //Vector Ta( n);
+        ldv_(Ta, Ld,a ,&n);
+//        Ta2 = tv_(T1,a.base(),&n);
+        tv_(TV, T1 ,a,&n);
+//        Ta.add(1.0,ConstVector(Ta2.base(), n));
+//        Ta.add(1.0,TV);
+          for (int j=0; j<n;++j)
+            Ta[j]+=TV[j];
+
+
+#ifdef TIMING_LOOP
+        }
+      mexPrintf("QT array mvm: finished: %d loops\n",loops);
+#endif
+      // create output and upload output data
+      /*      if (nlhs >= 1)
+      {
+      plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
+      double * mxinv= mxGetPr(plhs[0]);
+      // allocate likelihood array
+      for (int j=0;j<nper;++j)
+      mxinv[j]=(*vll)[j];
+      }
+      */   
+      mxFree(T1);
+      mxFree(Ld);
+      mxFree(TV);
+      } 
+    catch (const TSException& e) 
+      {
+      mexErrMsgTxt(e.getMessage());
+      } 
+    catch (SylvException& e) 
+      {
+      char mes[300];
+      e.printMessage(mes, 299);
+      mexErrMsgTxt(mes);
+      }
+    
+    } // mexFunction
+  }; // extern 'C'

mex/sources/kalman/matlab/qtmvm.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 QT f90 library by Andrea Pagano 
+% to multiply Quasi trinagular matrix (T) with a vector a
+%
+% function [a] = qtmvm(QT,a)
+%
+%  1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
+%  2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+%    a                      [double]    mm state vector.
+%
+% OUTPUTS
+%    Tinverse               [double]    mm*mm transition matrix of the state equation.
+**********************************************************/
+
+#include "qt.h"
+#include "kalman.h"
+#include "ts_exception.h"
+
+#include "GeneralMatrix.h"
+#include "Vector.h"
+#include "SylvException.h"
+
+#include "mex.h"
+
+extern "C" {
+  void mexFunction(int nlhs, mxArray* plhs[],
+    int nrhs, const mxArray* prhs[])
+    {
+    
+    if (nrhs < 2 )
+      mexErrMsgTxt("Must have  min 2 input parameters.\n");
+    if (nlhs < 1 )
+      mexErrMsgTxt("Must have  min 1 output parameters.\n")
+      
+      ;
+    //int start = 1; // default start of likelihood calculation
+    // test for univariate case
+    try 
+      {
+      // make input matrices
+      int n=mxGetM(prhs[0]);
+      
+      ConstGeneralMatrix QT(mxGetPr(prhs[0]), n, mxGetN(prhs[0]));
+      //      ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      Vector a (mxGetPr(prhs[1]), n);
+      
+      // create output and upload output data
+      plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
+      //      double * mxinv= mxGetPr(plhs[0]);
+      //      GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
+      Vector Ta (mxGetPr(plhs[0]), n);
+      //    Tinv.unit();
+      //      Ta.zeros();
+      
+        GeneralMatrix T1gm(n,n);
+        GeneralMatrix Ld(n,n);
+        Vector TV( n);
+      
+      // make storage for output
+      
+#ifdef TIMING_LOOP
+      int loops=1;//000;
+      if (nrhs >2 )
+        loops = (int)mxGetScalar(prhs[2]);
+      for (int tt=0;tt<loops;++tt)
+        {
+#endif
+#ifdef DEBUG	
+      //  QT.print();
+#endif
+        // 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
+     //   double *T1, *Ld, *dTa;//, dT1=-7.77;
+//  mexPrintf("start dT1 = %f\n", dT1);
+//        dT1 = qt2t_(QT.base() ,&n) ;
+// mexPrintf("end dT1 = %f\n", dT1);