SASification of sci_lu
Bernard HUGUENEY [Tue, 10 Mar 2009 08:30:45 +0000 (09:30 +0100)]
scilab/modules/linear_algebra/Makefile.am
scilab/modules/linear_algebra/includes/lu.h [new file with mode: 0644]
scilab/modules/linear_algebra/sci_gateway/c/sci_inv.c
scilab/modules/linear_algebra/sci_gateway/c/sci_lu.c
scilab/modules/linear_algebra/src/c/lu.c [new file with mode: 0644]

index 568154a..b12fb0d 100644 (file)
@@ -5,6 +5,7 @@
 
 LINEAR_ALGEBRA_C_SOURCES = src/c/schurtable.c \
 src/c/invert_matrix.c \
+src/c/lu.c \
 src/c/issymmetric.c
 
 LINEAR_ALGEBRA_FORTRAN_SOURCES = src/fortran/intdggbal.f \
diff --git a/scilab/modules/linear_algebra/includes/lu.h b/scilab/modules/linear_algebra/includes/lu.h
new file mode 100644 (file)
index 0000000..4f1f0ab
--- /dev/null
@@ -0,0 +1,33 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) ????-2009 - INRIA
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution.  The terms
+ * are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+#ifndef LU_H
+#define LU_H
+
+/**
+ * iLuM performs memory allocations and computes (calling iLu) LU decomposition 
+ * in :
+ * @param pData double* pointer to memory (either in 'd' or 'z' format) of argument matrix
+ * @param iRows int number of rows 
+ * @param iCols int number of columns 
+ * @param complexArg int (really bool) tells if data is in real (in 'd' format) or complex (in 'z' format)
+ * out :
+ * @param pdblLData double* pointer to memory (either in 'd' or 'z' format) where L result will be stored
+ * @param pdblUData double* pointer to memory (either in 'd' or 'z' format) where U result will be stored
+ * @param pdblEData double* NULL if E does not need to be computed, otherwise pointer to memory
+ * (either in 'd' or 'z' format) where E result will be stored
+ * @return if the operation successed (true) or not (false)
+ */
+extern int iLuM(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pbdlEData );
+extern int iLu(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pbdlEData
+              , int* piPivot, int* piWork, double* pdblWork);
+
+#endif
index d02a66d..03b1edd 100644 (file)
@@ -1,7 +1,7 @@
 
 /*
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) ????-2008 - INRIA
+ * Copyright (C) ????-2009 - INRIA
  *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
@@ -74,7 +74,7 @@ int C2F(intinv)(char *fname,unsigned long fname_len)
        if(ret >0){ Error(ret) ; }
        else{
          LhsVar(1) = 1;
-         // TODO rajouter le PutLhsVar(); quand il sera enlevé du gw_
+         /* TODO rajouter le PutLhsVar(); quand il sera enlevé du gw_ */
        }
        }
        return 0;
index e54c600..40f0a76 100644 (file)
@@ -1,7 +1,7 @@
 
 /*
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) ????-2008 - INRIA
+ * Copyright (C) ????-2009 - INRIA
  *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
  */
 
 #include <string.h>
-#include <stdio.h>
 #include "stack-c.h"
+#include "MALLOC.h"
 #include "gw_linear_algebra.h"
 #include "Scierror.h"
 #include "localization.h"
-/*--------------------------------------------------------------------------*/
-extern int C2F(intdgetrf)(char *fname, unsigned long fname_len);
-extern int C2F(intzgetrf)(char *fname, unsigned long fname_len);
-/*--------------------------------------------------------------------------*/
+#include "lu.h"
+
+#include "stack2.h" /* to create eye matrix */
+
 int C2F(intlu)(char *fname,unsigned long fname_len)
 {
-       int *header1;
-       int CmplxA;
-       int ret;
-
-       /*   lu(A)  */
-       if (GetType(1)!=sci_matrix) 
-       {
-               OverLoad(1);
-               return 0;
-       }
-       header1 = (int *) GetData(1);
-       CmplxA = header1[3];
-       switch (CmplxA) 
-       {
-               case REAL:
-                       ret = C2F(intdgetrf)("lu",2L);
-               break;
+  int complexArg;
+  int ret;
+  int iRows, iCols;
+  double* pData;
+  double* pDataReal;
+  double* pDataImg;
+  
+  /*   lu(A)  */
+  if ( (Rhs >=1 ) && GetType(1)!=sci_matrix) 
+    {
+      OverLoad(1);
+      return 0;
+    }
+  CheckRhs(1,1); /* one and only one arg */
+  CheckLhs(2,3); /* [L,U,[E]] = lu(A) */
+  complexArg=iIsComplex(1);
+  if(complexArg)
+    {
+      GetRhsVarMatrixComplex(1, &iRows, &iCols, &pDataReal, &pDataImg);
+      /* c -> z */
+      pData=(double*)oGetDoubleComplexFromPointer( pDataReal, pDataImg, iRows * iCols);
+    }
+       else
+         {
+           GetRhsVarMatrixDouble(1, &iRows, &iCols, &pData);
+         }
+       if( (iCols == 0) || (iRows == 0))
+         {
+           double* pdblL= NULL;
+           LhsVar(1)= 1;
+           ret = iAllocMatrixOfDouble(2, 0, 0, &pdblL);
+           /* check pbdlL or ret */
+           LhsVar(2)= 2;
+           if(Lhs == 3)
+             {
+               double* pdblE= NULL;
+               ret = iAllocMatrixOfDouble(2, 0, 0, &pdblE);
+               /* check pbdlE or ret */
+               LhsVar(3)= 3;
+             }
+         }
+       else
+         {
+           if( (iCols == -1) && (iRows == -1)) /* Rhs(1)=k*eye() => Lhs(1)=eye() Lhs(2)=k*eye(), Lhs(3)=eye() */
+             {
+               LhsVar(1)= 1;
+               if(complexArg)
+                 {
+                   double* pdblLReal;
+                   double* pdblLImg;
+                   iAllocComplexMatrixOfDouble(2, -1, -1, &pdblLReal, &pdblLImg);
+                   *pdblLReal= *pDataReal;
+                   *pdblLImg= *pDataImg;
+                   *pDataReal= 1.;
+                   *pDataImg= 0.;
+                 }
+               else
+                 {
+                   double* pdblLData;
+                   iAllocMatrixOfDouble(2, -1, -1, &pdblLData);
+                   *pdblLData= *pData;
+                   *pData= 1.;
+               }
+               LhsVar(2)= 2;
+               if(Lhs == 3)
+                 {
+                   if(complexArg)
+                     {
+                       double* pdblEReal;
+                       double* pdblEImg;
+                       iAllocComplexMatrixOfDouble(3, -1, -1, &pdblEReal, &pdblEImg);
+                       *pdblEReal= 1.;
+                       *pdblEImg= 0.;
+                     }
+                   else
+                     {
+                       double* pdblEData;
+                       iAllocMatrixOfDouble(3, -1, -1, &pdblEData);
+                       *pdblEData= 1.;
+                     }
+                   LhsVar(3)= 3;
+                 }
+             }
+           else
+             {
+               double *pdblLData;
+               double *pdblLReal;
+               double *pdblLImg;
+               double *pdblUData;
+               double *pdblUReal;
+               double *pdblUImg;
+               double *pdblEData;
+               int iMinRowsCols;
 
-               case COMPLEX:
-                       ret = C2F(intzgetrf)("lu",2L);
-               break;
+               pdblEData= NULL;
+               iMinRowsCols= Min(iRows, iCols);
 
-               default:
-                       Scierror(999,_("%s: Wrong type for input argument #%d: Real or Complex matrix expected.\n"),
-                       fname,1);
-               break;
-       }
+               if(complexArg)
+                 {
+                   iAllocComplexMatrixOfDouble(2, iRows, iMinRowsCols, &pdblLReal, &pdblLImg);
+                   iAllocComplexMatrixOfDouble(3, iMinRowsCols, iCols, &pdblUReal, &pdblUImg);
+                   /*
+                     we can allocate matrix of 'z' instead of calling oGetDoubleComplexFromPointer because the freshly allocated
+                     complex matrix does not contain any useful data.
+                    */
+                   pdblLData = (double*)MALLOC(iRows * iMinRowsCols * sizeof(doublecomplex) );
+                   pdblUData = (double*)MALLOC(iMinRowsCols * iCols * sizeof(doublecomplex) );
+                 }
+               else
+                 {
+                   iAllocMatrixOfDouble(2, iRows, iMinRowsCols, &pdblLData);
+                   iAllocMatrixOfDouble(3, iMinRowsCols, iCols, &pdblUData);
+                 }
+               if(Lhs == 3)
+                 {
+                   iAllocMatrixOfDouble(4, iRows, iRows, &pdblEData);
+                 }
+               ret =iLuM(pData, iRows, iCols, complexArg, pdblLData, pdblUData, pdblEData );
+               if(complexArg)
+                 {
+                   vGetPointerFromDoubleComplex((doublecomplex*)pdblLData, iRows * iMinRowsCols, pdblLReal, pdblLImg);
+                   FREE(pdblLData);
+                   vGetPointerFromDoubleComplex((doublecomplex*)pdblUData, iMinRowsCols * iCols, pdblUReal, pdblUImg);
+                   FREE(pdblUData);
+                 }
+               LhsVar(1)= 2;
+               LhsVar(2)= 3;
+               if(Lhs == 3)
+                 {
+                   LhsVar(3)= 4;
+                 }
+             }
+         }
+         /* TODO rajouter le PutLhsVar(); quand il sera enlevé du gw_ */
        return 0;
 }
 /*--------------------------------------------------------------------------*/
diff --git a/scilab/modules/linear_algebra/src/c/lu.c b/scilab/modules/linear_algebra/src/c/lu.c
new file mode 100644 (file)
index 0000000..0a64e6f
--- /dev/null
@@ -0,0 +1,168 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2009-2009 - DIGITEO - Bernard HUGUENEY
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+#include "machine.h"
+#include "core_math.h"
+#include "MALLOC.h"
+#include "lu.h"
+#include "doublecomplex.h"
+
+extern int F2C(zgetrf)(int*, int*, double*, int*, int*, int*);
+extern int F2C(dgetrf)(int*, int*, double*, int*, int*, int*);
+
+int iLuM(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pdblEData )
+{
+  int ret;
+  int* piPivot;
+  int* piWork;
+  int iMinRowsCols;
+  double* pdblWork;
+  iMinRowsCols= Min( iRows, iCols);
+  piPivot= (int*) MALLOC(iMinRowsCols * sizeof(int));
+  piWork= NULL;
+  pdblWork= NULL;
+  if(pdblEData == NULL) /* must allocate more tmp memory when not computing E */
+    {
+      piWork= (int*) MALLOC(iRows * sizeof(int));
+      pdblWork= (double*) MALLOC(iRows * iMinRowsCols * ( complexArg ? sizeof(doublecomplex) : sizeof(double)));
+    }
+  ret = iLu(pData, iRows, iCols, complexArg, pdblLData, pdblUData, pdblEData, piPivot, piWork, pdblWork);
+  if(pdblEData == NULL)
+    {
+      FREE( piWork);
+      FREE( pdblWork );
+    }
+  return ret;
+}
+int iLu(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pdblEData
+       , int* piPivot, int* piWork, double* pdblWork)
+{
+  int ret;
+  int iInfo;
+  if(complexArg)
+    {
+      C2F(zgetrf)(&iRows, &iCols, pData, &iRows, piPivot, &iInfo);
+    }
+  else
+    {
+      C2F(dgetrf)(&iRows, &iCols, pData, &iRows, piPivot, &iInfo);
+    }
+  if (iInfo <0)
+    {
+      ret= iInfo;
+    }
+  else
+    {
+      int iMinRowsCols;
+      int i, j;
+      double* pdCurrentL;
+      double* pdCurrentU;
+      double* pdCurrentArg;
+      iMinRowsCols = Min(iRows, iCols);
+      /* fill L matrix with 1.[+0i] on the diagonal */
+      /*
+       /!\ if E was not requested, we will have to swap rows before returning the result, so we do not
+       fill the real L matrix, but pdblWork instead. We will then copy and swap rows from pdblWork to the result matrix.
+       */
+      pdCurrentL= (pdblWork == NULL) ? pdblLData : pdblWork ;
+      pdCurrentArg= pData;
+      for( j= 0; j != iMinRowsCols; ++j)
+       {
+         for ( i=0; i!= iRows; ++i, ++pdCurrentL, ++pdCurrentArg)
+           {
+             if(i > j)
+               {
+                 *pdCurrentL = *pdCurrentArg;
+                 if(complexArg)
+                   { /* copy imaginary part */
+                     *++pdCurrentL = *++pdCurrentArg;
+                   }
+               }
+             else
+               {
+                 *pdCurrentL= (i==j) ? 1. : 0.;
+                 if(complexArg)
+                   { /* set imaginary part to 0. and advance pdCurrentArg */
+                     *++pdCurrentL= 0.;
+                     ++pdCurrentArg;
+                   }
+               }
+           }
+       }
+      /* init U */
+      pdCurrentU= pdblUData;
+      pdCurrentArg= pData;
+      for( j= 0; j != iCols; ++j)
+       {
+         pdCurrentArg= pData + j * iMinRowsCols * (complexArg ? 2 : 1);
+         for ( i=0; i!= iMinRowsCols; ++i, ++pdCurrentU, ++pdCurrentArg)
+           {
+             if( i <= j)
+               {
+                 *pdCurrentU = *pdCurrentArg;
+                 if(complexArg)
+                   { /* copy imaginary part */
+                     *++pdCurrentU = *++pdCurrentArg;
+                   }
+               }
+             else
+               {
+                 *pdCurrentU= 0.;
+                 if(complexArg)
+                   {
+                     *++pdCurrentU= 0.;
+                     ++pdCurrentArg;
+                   }
+               }
+           }
+       }
+      if( pdblEData == NULL)
+       {
+         for(i=0; i != iRows; ++i)
+           {
+             piWork[i]=i;
+           }
+         for(i=0; i != iMinRowsCols; ++i)
+           {
+             int ip;
+             ip= piPivot[i]-1; /* -1 because piPivot contains Fortran index*/
+             if( ip != i) 
+               {
+                 int swapTmp;
+                 swapTmp= piWork[i];
+                 piWork[i]= piWork[ip];
+                 piWork[ip]= swapTmp;
+               }
+           }
+         for(i= 0; i != iRows; ++i)
+           {
+             int ip;
+             ip=piWork[i];
+             if(complexArg)
+               {
+                 C2F(zcopy)(&iCols, pdblWork+(i*2), &iRows, pdblLData+(ip*2), &iRows); 
+               }
+             else
+               {
+                 C2F(dcopy)(&iCols, pdblWork+i, &iRows, pdblLData+ip, &iRows);
+               }
+           }
+       }/* end if E was not requested */
+      else
+       {
+         double dblZero=0., dblOne=1.;
+         int iOne=1;
+         F2C(dlaset)("F", &iRows, &iRows, &dblZero, &dblOne, pdblEData, &iRows);
+         F2C(dlaswp)(&iRows, pdblEData, &iRows, &iOne, &iMinRowsCols, piPivot, &iOne);
+       }
+    }/* end if iInfo reported an error in [z|d]getrf */
+  return ret;
+}