SASification of linear_algebra sci_eigen.c for Scilab spec() function
Bernard HUGUENEY [Wed, 25 Mar 2009 08:28:25 +0000 (09:28 +0100)]
scilab/modules/linear_algebra/Makefile.am
scilab/modules/linear_algebra/Makefile.in
scilab/modules/linear_algebra/includes/eigen.h [new file with mode: 0644]
scilab/modules/linear_algebra/sci_gateway/c/sci_eig.c
scilab/modules/linear_algebra/src/c/eigen.c [new file with mode: 0644]

index 7368ef5..8b60388 100644 (file)
@@ -8,6 +8,7 @@ src/c/invert_matrix.c \
 src/c/lu.c \
 src/c/qr.c \
 src/c/hess.c \
+src/c/eigen.c \
 src/c/issymmetric.c
 
 LINEAR_ALGEBRA_FORTRAN_SOURCES = src/fortran/intdggbal.f \
index 5eb416a..46414c5 100644 (file)
@@ -113,6 +113,7 @@ am__objects_2 = libscilinear_algebra_la-schurtable.lo \
        libscilinear_algebra_la-invert_matrix.lo \
        libscilinear_algebra_la-lu.lo libscilinear_algebra_la-qr.lo \
        libscilinear_algebra_la-hess.lo \
+       libscilinear_algebra_la-eigen.lo \
        libscilinear_algebra_la-issymmetric.lo
 am__objects_3 = SELECT.lo Ex-schur.lo
 am__objects_4 = libscilinear_algebra_la-sci_backslash.lo \
@@ -378,6 +379,7 @@ src/c/invert_matrix.c \
 src/c/lu.c \
 src/c/qr.c \
 src/c/hess.c \
+src/c/eigen.c \
 src/c/issymmetric.c
 
 LINEAR_ALGEBRA_FORTRAN_SOURCES = src/fortran/intdggbal.f \
@@ -626,6 +628,7 @@ distclean-compile:
        -rm -f *.tab.c
 
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscilinear_algebra_la-assembleEigenvectors.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscilinear_algebra_la-eigen.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscilinear_algebra_la-gw_linear_algebra.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscilinear_algebra_la-gw_linear_algebra2.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscilinear_algebra_la-hess.Plo@am__quote@
@@ -712,6 +715,13 @@ libscilinear_algebra_la-hess.lo: src/c/hess.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscilinear_algebra_la_CFLAGS) $(CFLAGS) -c -o libscilinear_algebra_la-hess.lo `test -f 'src/c/hess.c' || echo '$(srcdir)/'`src/c/hess.c
 
+libscilinear_algebra_la-eigen.lo: src/c/eigen.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscilinear_algebra_la_CFLAGS) $(CFLAGS) -MT libscilinear_algebra_la-eigen.lo -MD -MP -MF $(DEPDIR)/libscilinear_algebra_la-eigen.Tpo -c -o libscilinear_algebra_la-eigen.lo `test -f 'src/c/eigen.c' || echo '$(srcdir)/'`src/c/eigen.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscilinear_algebra_la-eigen.Tpo $(DEPDIR)/libscilinear_algebra_la-eigen.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/eigen.c' object='libscilinear_algebra_la-eigen.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscilinear_algebra_la_CFLAGS) $(CFLAGS) -c -o libscilinear_algebra_la-eigen.lo `test -f 'src/c/eigen.c' || echo '$(srcdir)/'`src/c/eigen.c
+
 libscilinear_algebra_la-issymmetric.lo: src/c/issymmetric.c
 @am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscilinear_algebra_la_CFLAGS) $(CFLAGS) -MT libscilinear_algebra_la-issymmetric.lo -MD -MP -MF $(DEPDIR)/libscilinear_algebra_la-issymmetric.Tpo -c -o libscilinear_algebra_la-issymmetric.lo `test -f 'src/c/issymmetric.c' || echo '$(srcdir)/'`src/c/issymmetric.c
 @am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscilinear_algebra_la-issymmetric.Tpo $(DEPDIR)/libscilinear_algebra_la-issymmetric.Plo
diff --git a/scilab/modules/linear_algebra/includes/eigen.h b/scilab/modules/linear_algebra/includes/eigen.h
new file mode 100644 (file)
index 0000000..87c0d76
--- /dev/null
@@ -0,0 +1,110 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) ????-2009 - INRIA 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
+ *
+ */
+#ifndef EIGEN_H
+#define EIGEN_H
+
+#include "doublecomplex.h"
+
+/*
+ * Functions used to compute eigen values (and possibly eigen vectors) of one [complex|real] [un]symmetric matrix.
+ * We use 4 different functions because for e;fficiency reasons, API had to be different :
+ * 1) eigen vectors for symmetric matrix are computed in place in the same memory location as the argument matrix
+ * 2) complex eigen values for unsymmetric matrix are computed either in 'c' format (for real argument matrix) or 'z' format (complex argument matrix)
+ *
+ * @param pData double[complex]* in (& out for symmetric) argument matrix data, for symmetric, is also use to output eigen vectors
+ * @param iCols int nb of rows/cols of matrix
+ * @param computeEigenVectors int (boolean semantics) only need for symmetric : indicates weither eigenvectors should be computed
+ *
+ * out :
+ * eigen values either real or complex in 'z' format
+ * @param pEigenValues double[complex]* output eigen values
+ *
+ *  or complex in 'c' format
+ * @param pEigenValuesReal double* output eigen values real part
+ * @param pEigenValuesImg double* output eigen values imaginary part
+ *
+ * for unsymmetric cases, output complex eigen vectors either in 'z' format (complex argument)
+ * @param pEigenVectors doublecomplex* output eigen vectors
+ *
+ * or in 'c' format (real argument)
+ * @param pEigenVectorsReal double* output eigen vectors real part
+ * @param pEigenVectorsImg double* output eigen vectors imaginary part
+ *
+ */
+extern int iEigen1ComplexSymmetricM(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues);
+extern int iEigen1RealSymmetricM(double* pData, int iCols, int computeEigenVectors, double* pEigenValues);
+
+extern int iEigen1ComplexM(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pEigenVectors);
+extern int iEigen1RealM(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pEigenVectorsReal, double* pEigenVectorsImg);
+
+/*
+ * With 2 lhs, eigenvalues are expanded from a 1D vector to the diagonal of a matrix.
+ * The resulting matrix is a diagonal one: m[i,j]=0. for i!=j .
+ * For efficiency reasons, we also convert from 'z' to 'c' format at the same time. In this case, conversion cannont be done in place.
+ *
+ * /!\ For the real case, expansion is done in place, so obviously the memory must have been reserved accordingly /!\
+ * (i.e for the iCols x iCols matrix, not only the iCols x 1 vector.)
+ *
+ * @param pData double[complex]* : input  (&output for real) pData[0:iCols[ contains diagonal elements
+ * @param iCols int : input nb of elements
+ * 
+ * out:
+ * @param  pRMatrix double* real part of the diagonal output matrix
+ * @param  pIMatrix double* imaginary part of the diagonal output matrix
+ */
+extern void expandToDiagonalOfMatrix(double* pData, int iCols);
+extern void expandZToDiagonalOfCMatrix(doublecomplex const* pZVector, int iCols, double* pRMatrix, double* pIMatrix);
+
+
+
+/*
+ * Functions used to compute eigen values (and possibly eigen vectors) of two [complex|real]  matrix.
+ *
+ * @param pData1 double[complex]* in argument1 matrix data
+ * @param pData2 double[complex]* in argument2 matrix data
+ *
+ * @param iCols int nb of rows/cols of matrix
+ *
+ * out : when something should not be computed (depending on lhs), output buffer ptr is NULL
+ *
+ * Alpha is either
+ * in 'z' format for complex matrix
+ * @param pAlpha doublecomplex*
+ *
+ * or 'c' format for real matrix
+ *
+ * @param pAlphaReal double*
+ * @param pAlphaImg double*
+ *
+ *
+ * @param pBeta double[complex]* 
+ *
+ * L and R are either
+ * in 'z' format for complex matrix
+ * @param pL doublecomplex*
+ * @param pR doublecomplex*
+ *
+ * or 'c' format for real matrix
+ *
+ * @param pLReal double*
+ * @param pLImg double*
+ * @param pRReal double*
+ * @param pRImg double*
+ *
+ */
+
+extern int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols
+                    , doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL);
+extern int iEigen2RealM(double* pData1, double* pData2, int iCols
+                    , double* pAlphaReal, double* pAlphaImg, double* pBeta, double* pRReal, double* pRimg, double* pLReal, double* pLImg);
+
+#endif
index 6f274a7..7425c73 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 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
 
 
 #include <string.h>
-#include <stdio.h>
 #include "stack-c.h"
 #include "issymmetric.h"
 #include "gw_linear_algebra.h"
 #include "Scierror.h"
 #include "localization.h"
-/*--------------------------------------------------------------------------*/
+#include "vfinite.h"
+#include "MALLOC.h"
+#include "sciprint.h"
+#include "eigen.h"
+/*
+evals=spec(A)
+[R,diagevals]=spec(A)
+
+evals=spec(A,B)
+[alpha,beta]=spec(A,B)
+[alpha,beta,Z]=spec(A,B)
+[alpha,beta,Q,Z]=spec(A,B)
+*/
 extern int C2F(complexify)(int *num);
-/*--------------------------------------------------------------------------*/
-int C2F(inteig)(char *fname,unsigned long fname_len)
+
+/* internal function used to detect imaginary part completly and strictly equal to 0. causing troubles to LAPACK. */  
+static int isArrayZero(int size , double * values)
 {
-  int *header1, *header2;
-  int CmplxA, CmplxB;
-  int ret;
-  int Symmetric;
-  int X;
+  while(--size)
+    {
+      if(*values!=0.) /* not a bug: only 0. is of interest, not epsilon */
+       {
+         return 0;
+       }
+    }
+  return 1;
+}
 
-  switch (Rhs) 
-  {
-       case 1:   /* spec(A)   */
-               if (GetType(1)!=sci_matrix) 
+int C2F(inteig)(char *fname,unsigned long fname_len)
+{
+  int ret=0;
+  CheckRhs(1,2);
+  CheckLhs(1,4);
+  if(GetType(1) != sci_matrix)
+    {
+      OverLoad(1);
+      ret= 0;
+    }
+  else
+    {
+      switch (Rhs) {
+      case 1:{
+       int iRows, iCols;
+       doublecomplex* pData= NULL;
+       double* pDataReal= NULL;
+       double* pDataImg= NULL;
+       int complexArg= iIsComplex(1);
+       CheckLhs(1,2);
+       if(complexArg)
+         {
+           GetRhsVarMatrixComplex(1, &iRows, &iCols, &pDataReal, &pDataImg);
+           /* c -> z */
+           pData=oGetDoubleComplexFromPointer( pDataReal, pDataImg, iRows * iCols);
+           if(!pData)
+             {
+               Scierror(999,_("%s: Cannot allocate more memory.\n"),fname);
+               ret = 1;
+             }
+         }
+       else
+         {
+           GetRhsVarMatrixDouble(1, &iRows, &iCols, &pDataReal);
+         }
+       if( iCols == -1)
+         {
+           Scierror(999,_("Size varying argument a*eye(), (arg %d) not allowed here.\n"), 1);
+         }
+       else /* not eye() matrix */
+         {
+           if (iCols==0)
+             {
+               if(Lhs == 1)
+                 {
+                   LhsVar(1)= 1;
+                 }
+               else /* Lhs == 2 */
+                 {
+                   double* dummy;
+                   int unused = complexArg 
+                     ? iAllocComplexMatrixOfDouble(2, iCols, iCols, &dummy, &dummy)
+                     : iAllocMatrixOfDouble(2, iCols, iCols, &dummy);
+                   LhsVar(1)= 2;
+                   LhsVar(2)= 1;
+                 }
+             }
+           else
+             {/* now the interesting case, can be either complex or real and symmetric or not */
+           
+               double* pEigenValuesReal= NULL;
+               double* pEigenValuesImg= NULL;
+               double* pEigenVectorsReal= NULL;
+               double* pEigenVectorsImg= NULL;
+               int symmetric= 0; /* bool in fact */
+               int const eigenValuesCols= (Lhs==1) ? 1 : iCols ;
+               int const totalSize= iRows * iCols;
+               if ( !(complexArg 
+                      ? C2F(vfiniteComplex)(&totalSize, pData) 
+                      : C2F(vfinite)(&totalSize, pDataReal)))
+                 {
+                   SciError(264);
+                   return 0;
+                 }
+               if( (symmetric=C2F(issymmetric)(&Rhs) ) )
+                 {
+                   iAllocMatrixOfDouble(2, iCols, eigenValuesCols, &pEigenValuesReal);
+                   /* if matrix is symmetric, the eigenvectors can reuse Rhs because the matrix is of the same type & dimensions */
+                 }
+               else
+                 {
+                   iAllocComplexMatrixOfDouble(2, iCols, eigenValuesCols, &pEigenValuesReal, &pEigenValuesImg);
+                   if(Lhs==2)
+                     {
+                       iAllocComplexMatrixOfDouble(3, iCols, iCols, &pEigenVectorsReal, &pEigenVectorsImg);
+                     }
+                 }
+               if( complexArg )
+                 {
+                   if(symmetric)
+                     {
+                       ret= iEigen1ComplexSymmetricM(pData, iCols, (Lhs==2), pEigenValuesReal);
+                       if(Lhs == 2)
+                         { /* we reuse memory from rhs1 to put back the resulting eigenvectors */
+                           vGetPointerFromDoubleComplex(pData, iCols * iCols , pDataReal, pDataImg);
+                           expandToDiagonalOfMatrix(pEigenValuesReal, iCols);
+                           LhsVar(1)= 1;
+                           LhsVar(2)= 2;
+                         }
+                       else
+                         {
+                           LhsVar(1)= 2;
+                         }
+                     }
+                   else
+                     {
+                       doublecomplex* pEigenValues= (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex));
+                       doublecomplex* pEigenVectors = (Lhs== 2) ? (doublecomplex*)MALLOC(sizeof(doublecomplex) * iCols * iCols) : NULL ;
+                       ret= iEigen1ComplexM(pData, iCols, pEigenValues, pEigenVectors);
+                       if(Lhs==2)
+                         {
+                           expandZToDiagonalOfCMatrix(pEigenValues, iCols, pEigenValuesReal, pEigenValuesImg); 
+                           vGetPointerFromDoubleComplex(pEigenVectors, iCols * iCols, pEigenVectorsReal, pEigenVectorsImg);
+                           FREE(pEigenVectors);
+                           LhsVar(1)= 3;
+                           LhsVar(2)= 2;
+                         }
+                       else
+                         {
+                           vGetPointerFromDoubleComplex(pEigenValues, iCols, pEigenValuesReal, pEigenValuesImg);
+                           LhsVar(1)= 2;
+                         }
+                       FREE(pEigenValues);
+                     }
+                 }
+               else /* real */
+                 {
+                   if(symmetric)
+                     {
+                       ret= iEigen1RealSymmetricM(pDataReal, iCols, (Lhs==2), pEigenValuesReal);
+                       if(Lhs == 2)
+                         {
+                           expandToDiagonalOfMatrix(pEigenValuesReal, iCols);
+                           LhsVar(1)= 1;
+                           LhsVar(2)=2;
+                         }
+                       else
+                         {
+                           LhsVar(1)=2;
+                         }
+                     }
+                   else
+                     {
+                       ret= iEigen1RealM(pDataReal, iCols, pEigenValuesReal, pEigenValuesImg, pEigenVectorsReal, pEigenVectorsImg);
+                       if( Lhs == 2)
+                         {
+                           expandToDiagonalOfMatrix(pEigenValuesReal, iCols);
+                           expandToDiagonalOfMatrix(pEigenValuesImg, iCols);
+                           LhsVar(1)= 3;
+                           LhsVar(2)= 2;
+                         }
+                       else
+                         {
+                           LhsVar(1)= 2;
+                         }
+                     }
+                 }
+             }
+         }
+       if(complexArg)
+         {
+           vFreeDoubleComplexFromPointer(pData);
+         }
+       break;
+      }
+      case 2:{
+       if(GetType(2) != sci_matrix)
+         {
+           OverLoad(2);
+           ret= 0;
+         }
+       {
+         doublecomplex* pData[2];
+         double* pDataReal[2];
+         double* pDataImg[2];
+         int iCols[2], iRows[2];
+         int i;
+         /* Because of bug #3652 ( http://bugzilla.scilab.org/show_bug.cgi?id=3652 ), complex matrix with both imaginary part null
+          * must be handled as real matrix.
+          */
+         int complexArgs= iIsComplex(1) || iIsComplex(2);
+         /* if either matrix is complex, promote the other to complex */
+         if(complexArgs)
+           {
+             int byRef=1;
+             C2F(complexify)(&byRef);
+             byRef=2;
+             C2F(complexify)(&byRef);
+           }
+         if(complexArgs)
+           {
+             for(i=0; i!=2; ++i)
                {
-                       OverLoad(1);
-                       return 0;
+                 GetRhsVarMatrixComplex(1+i, &iRows[i], &iCols[i], &pDataReal[i], &pDataImg[i]);
+                 if( !(pData[i]= oGetDoubleComplexFromPointer( pDataReal[i], pDataImg[i], iRows[i] * iCols[i])) )
+                   {
+                     Scierror(999,_("%s: Cannot allocate more memory.\n"),fname);
+                     ret = 1;
+                   }
                }
-               header1 = (int *) GetData(1);
-               CmplxA=header1[3];
-               X = 1;
-               Symmetric = C2F(issymmetric)(&X);
-               switch (CmplxA) 
+           }
+         else
+           {
+             for(i=0; i!=2; ++i)
                {
-                       case REAL:
-                               switch (Symmetric) 
-                               {
-                                       case NO :
-                                               ret = sci_dgeev("spec",4L);
-                                       break;
-                                       case YES :
-                                               ret = sci_dsyev("spec",4L);
-                                       break;
-                               }
-                       break;
-
-                       case COMPLEX:
-                               switch (Symmetric) 
-                               {
-                                       case NO :
-                                               ret = sci_zgeev("spec",4L);
-                                       break;
-                                       case YES:
-                                               ret = sci_zheev("spec",4L);
-                                       break;
-                               }
-                       break;
-
-                       default:
-                               Scierror(999,_("%s: Wrong type for input argument #%d: Real or Complex matrix expected.\n"),
-                                       fname,1);
-                       break;
-               } /* end switch  (CmplxA) */
-    break; /* end case 1 */
-
-       case 2: /* gspec(A,B) */
-               if (GetType(1)!=sci_matrix) 
+                 GetRhsVarMatrixDouble(1+i, &iRows[i], &iCols[i], &pDataReal[i]);
+               }
+           }
+         for( i=0; i!=2; ++i)
+           {
+             if (iRows[i]!= iCols[i])
                {
-                       OverLoad(1);
-                       return 0;
+                 Scierror(999,_("%s: Wrong size for input argument #%d: A square matrix expected.\n"), fname, i+1);
                }
-               if (GetType(2)!=sci_matrix) 
+           }
+         if( iCols[0] != iCols[1])
+           {/* /!\ reusing existing error msg, but it could be more explicit :"%s: arguments %d and %d must have equal dimensions.\n" */
+             Scierror(999,_("%s and %s must have equal dimensions.\n"),"A","B");
+           }
+         if( iCols[0] == 0)
+           {
+             switch (Lhs){
+             case 4:{
+               double* dummy;
+               if(complexArgs)
+                 {
+                   iAllocComplexMatrixOfDouble(4, 0, 0, &dummy, &dummy);
+                 }
+               else
+                 {
+                   iAllocMatrixOfDouble(4, 0, 0, &dummy) ;
+                 }
+               LhsVar(4)= 4; /* no break */
+             }
+             case 3:{
+               double* dummy;
+               if(complexArgs)
+                 {
+                   iAllocComplexMatrixOfDouble(3, 0, 0, &dummy, &dummy);
+                 }
+               else
+                 {
+                   iAllocMatrixOfDouble(3, 0, 0, &dummy) ;
+                 }
+               LhsVar(3)= 3;
+             }/* no break */
+             case 2:
+               LhsVar(2)= 2; /* no break */
+             default: /* case 1: as we have done CheckLhs(1,4)*/
+               LhsVar(1)= 1; 
+             }
+           }
+         else
+           {
+             int inf=0;
+             int const totalSize= iCols[0] * iCols[0];
+             for(i=0; i!=2 && !inf; ++i)
                {
-                       OverLoad(2);
-                       return 0;
+                 if( (inf= (inf || !(complexArgs 
+                                     ? C2F(vfiniteComplex)(&totalSize, pData[i]) 
+                                     : C2F(vfinite)(&totalSize, pDataReal[i]) ) )))
+                   {/* /!\ reusing error msg, but could be more explicit ny being prefixed by %s for fname */
+                     Scierror(999,_("Wrong value for argument %d: Must not contain NaN or Inf.\n"),i+1);
+                   }
                }
-               header1 = (int *) GetData(1);
-               header2 = (int *) GetData(2);
-               CmplxA=header1[3];
-               CmplxB=header2[3];
-               switch (CmplxA) 
+             if(!inf)
                {
-                       case REAL:
-                               switch (CmplxB) 
-                               {
-                                       case REAL :
-                                       /* A real, Breal */
-                                               ret = sci_dggev("gspec",5L);
-                                       break;
-                                       case COMPLEX :
-                                       /* A real, B complex : complexify A */
-                                               X=1;
-                                               C2F(complexify)(&X);
-                                               ret = sci_zggev("gspec",5L);
-                                       break;
-                                       default:
-                                               Scierror(999,_("%s: Wrong type for input argument #%d: Real or Complex matrix expected.\n"),
-                                                       fname,2);
-                                       break;
-                               }
-                       break;
-                       case COMPLEX :
-                               switch (CmplxB) 
-                               {
-                                       case REAL :
-                                       /* A complex, B real : complexify B */
-                                               X=2;
-                                               C2F(complexify)(&X);
-                                               ret = sci_zggev("gspec",5L);
-                                       break;
-                                       case COMPLEX :
-                                       /* A complex, B complex */
-                                               ret = sci_zggev("gspec",5L);
-                                       break;
-                                       default:
-                                                       Scierror(999,_("%s: Wrong type for input argument #%d: Real or Complex matrix expected.\n"),
-                                                       fname,2);
-                                       break;
-                               }
-                       break;
-                       default :
-                               Scierror(999,_("%s: Wrong type for input argument #%d: Real or Complex matrix expected.\n"),
-                                                       fname,1);
-                               break;
-               } /*end  switch (CmplxA) */
-    break;/* end case 2 */
+                 doublecomplex* pL= NULL;
+                 double* pLReal= NULL;
+                 double* pLImg= NULL;
+                 doublecomplex* pR= NULL;
+                 double* pRReal= NULL;
+                 double* pRImg= NULL;
+                 doublecomplex* pAlpha= NULL;
+                 double* pAlphaReal= NULL;
+                 double* pAlphaImg= NULL;
+                 doublecomplex* pBeta= NULL;
+                 double* pBetaReal= NULL;
+                 double* pBetaImg= NULL;
+                 /*
+                   special case for complex matrix with null imaginary parts.
+                 */
+                 int complexArgsWithZeroImg= isArrayZero(totalSize, pDataImg[0]) &&  isArrayZero(totalSize, pDataImg[1]);
+                 /* API could be better wrt c or z format of L and R but with less code ruse from assembleEigenVectors */
+                 switch(Lhs){/* for Complex case, it could be possible to reuse var from Rhs1 & Rhs2 for Lhs 3 & Lhs4 */
+                 case 4:
+                   iAllocComplexMatrixOfDouble(6, iCols[0], iCols[0], &pLReal, &pLImg);  /* nobreak */
+                 case 3:
+                   iAllocComplexMatrixOfDouble(5, iCols[0], iCols[0], &pRReal, &pRImg); /* nobreak */
+                 case 2:{
+                   if(complexArgs && !complexArgsWithZeroImg)
+                     {
+                       pBeta= iAllocComplexMatrixOfDouble(4, iCols[0], 1, &pBetaReal, &pBetaImg)
+                         ? NULL /* var alloc failed, no need to try to alloc 'z' vector */
+                         : (doublecomplex*)MALLOC(iCols[0] * sizeof(doublecomplex));
+
+                     }
+                   else
+                     {
+                       iAllocMatrixOfDouble(4, iCols[0], 1, &pBetaReal);
+                     }
+                 } /* nobreak */
+                 default : /*case 1 as we have done CheckLhs(1,4) */
+                   { /* note that real matrix use alpha in 'c' format while complex matrix use it in 'z' format, cf dggev and zggev */
+                     pAlpha = ( iAllocComplexMatrixOfDouble(3, iCols[0], 1, &pAlphaReal, &pAlphaImg) 
+                                || (!complexArgs || complexArgsWithZeroImg) )
+                       ? NULL /* var alloc failed or we are handling real matrix */
+                       : (doublecomplex*) MALLOC(iCols[0] * sizeof(doublecomplex)) ;
+                   }
+                   /* nobreak */
+                 }
+                 ret= (complexArgs && !complexArgsWithZeroImg)
+                   ? iEigen2ComplexM(pData[0], pData[1], iCols[0], pAlpha, pBeta, pR, pL)
+                   : iEigen2RealM(pDataReal[0], pDataReal[1], iCols[0], pAlphaReal, pAlphaImg, pBetaReal, pRReal, pRImg, pLReal, pLImg) ;
+                 if(ret >0 )
+                   {
+                     sciprint(_("Warning :\n"));
+                     sciprint(_("Non convergence in the QZ algorithm.\n"));
+                     sciprint(_("The top %d  x %d blocks may not be in generalized Schur form.\n"), ret);
+                   }
+                 if( ret >= 0)
+                   {
+                     switch (Lhs) 
+                       {
+                       case 4:
+                         if(complexArgs && !complexArgsWithZeroImg)
+                           {
+                             vGetPointerFromDoubleComplex(pL, totalSize, pLReal, pLImg);
+                           }/* nobreak */
+                       case 3:
+                         {
+                           if(complexArgs && !complexArgsWithZeroImg)
+                             {
+                               vGetPointerFromDoubleComplex(pR, totalSize, pRReal, pRImg);
+                             }
+                           if(Lhs == 4)
+                             {
+                               LhsVar(4)= 5;
+                               LhsVar(3)= 6;
+                             }
+                           else /* Lhs == 3 */
+                             {
+                               LhsVar(3)= 5;
+                             }
+                         }/* nobreak */
+                       case 2:
+                         {
+                           if(complexArgs && !complexArgsWithZeroImg)
+                             {
+                               vGetPointerFromDoubleComplex(pBeta, iCols[0], pBetaReal, pBetaImg);
+                             }
+                           LhsVar(2)= 4; /* nobreak */
+                         }
+                       default: /* case 1 */
+                         {
+                           if(complexArgs)
+                             {
+                               vGetPointerFromDoubleComplex(pAlpha, iCols[0], pAlphaReal, pAlphaImg);
+                             }
+                           LhsVar(1)= 3;
+                         }
+                       }
+                         
+                   }
+                 /* adjust*/
+                 FREE(pAlpha);
+                 FREE(pBeta);
+               }
+             if(complexArgs)
+               {
+                 vFreeDoubleComplexFromPointer(pData[0]);
+                 vFreeDoubleComplexFromPointer(pData[1]);
+               }
 
-  }/* end switch (Rhs) */
+           }
+       }
+       
+      }
+       break;/* end case 2 */
+       
+      }/* end switch (Rhs) */
+    }
   return 0;
 }
 /*--------------------------------------------------------------------------*/
diff --git a/scilab/modules/linear_algebra/src/c/eigen.c b/scilab/modules/linear_algebra/src/c/eigen.c
new file mode 100644 (file)
index 0000000..a617a5c
--- /dev/null
@@ -0,0 +1,648 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) ????-2009 - INRIA 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 <string.h> /* memset */
+#include "core_math.h"
+#include "machine.h"
+#include "MALLOC.h"
+
+#include "eigen.h"
+
+/*
+ * Auxiliary functions to expand Lapack compact representation of eigen vectors into scilab matrix
+ * see implementation at the end of the file for API documentation
+ */
+static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg, 
+                                             double * EVRealSource, double * EVRealTarget, double * EVImgTarget);
+static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg);
+
+
+static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
+                  , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork);
+
+static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
+               , double* pR, double* pL, double* pWork, int iWorkSize);
+
+static double* allocateDggevWorkspace(int iCols, int* pWorksize);
+
+
+/*
+ * LAPACK Fortran functions used to do the real work for eigen values / vectors of one matrix.
+ */
+
+extern void C2F(zheev)(char const JOBZ[1]/*'N'|'V'*/, char const UPLO[1]/*'U'|'L'*/, int const* N, doublecomplex* A
+                       , int const* LDA, double* W, doublecomplex* WORK, int const* LWORK, double* RWORK, int* INFO);
+
+extern void C2F(zgeev)(char const JOBVL[1]/*'N'|'V'*/,char const JOBVR[1]/*'N'|'V'*/, int const* N, doublecomplex* A
+                       , int const* LDA, doublecomplex* W, doublecomplex* VL, int const* LDVL, doublecomplex* VR, int const* LDVR
+                       ,doublecomplex* WORK, int const* LWORK, double* RWORK, int* INFO);
+
+extern void C2F(dsyev)( char const JOBZ[1]/*'N'|'V'*/, char const UPLO[1]/*'U'|'L'*/, int const* N, double* A
+                        , int const* LDA, double* W, double* WORK, int const* LWORK, int* INFO);
+
+extern void C2F(dgeev)(char const JOBVL[1]/*'N'|'V'*/,char const JOBVR[1]/*'N'|'V'*/, int const* N, double* A
+                       , int const* LDA, double* WR, double* WI, double* VL, int const* LDVL, double* VR, int const* LDVR
+                       ,double* WORK, int const* LWORK, int* INFO);
+
+
+/*
+ * LAPACK Fortran functions used to do the real work for eigen values / vectors of two matrix.
+ */
+
+
+extern void C2F(dggev)(char const jobVL[1]/*'N'|'V'*/,char const jobVR[1]/*'N'|'V'*/, int const* n, double* a, int const* ldA
+                 ,double* b, int const* ldB, double* alphaR, double* alphaI, double* beta, double* vl, int const* ldVl, double* vr, int const* ldVr
+                 ,double* WORK, int const* LWORK, int* INFO);
+extern void C2F(zggev)(char const jobVL[1]/*'N'|'V'*/,char const jobVR[1]/*'N'|'V'*/, int const* n, doublecomplex* a, int const* ldA
+                 ,doublecomplex* b, int const* ldB, doublecomplex* alpha, doublecomplex* beta, doublecomplex* vl, int const* ldVl
+                 , doublecomplex* vr, int const* ldVr, doublecomplex* WORK, int const* LWORK, double * RWORK, int* INFO);
+
+
+/*
+ * external function to perform division of complex vectors
+ */
+extern int C2F(wwrdiv)(double *ar, double *ai, int const *ia, double *br, double *bi, int const *ib
+                      , double *rr, double *ri, int const *ir, int const *n, int *ierr);
+
+
+/*
+ * Wrappers querying optimal worsksize and computing minimal worksize for LAPACK functions
+ * @param iCols int nb of rows/columns of the matrix
+ * @param lhs int nb of expected outputs from the spec scilab function, only when useful to compute optimal worksize
+ * out:
+ *
+ * @param optWorkSize int* points to the optimal worksize
+ * @param minWorkSize int* points to the minimal worksize
+ */
+static int zheevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
+{
+ int info=0, query=-1;
+  doublecomplex opt;
+  C2F(zheev)("N", "U", &iCols, NULL, &iCols, NULL, &opt, &query,NULL, &info );
+  *optWorkSize= (int) opt.r;
+  *minWorkSize= Max(1, 2*iCols-1);
+  return info;
+}
+
+static int zgeevWorkSizes(int iCols,  int lhs, int* optWorkSize, int* minWorkSize)
+{
+  int info=0, query=-1;
+  doublecomplex opt;
+  C2F(zgeev)("N", (lhs == 1 ? "N" : "V"), &iCols, NULL, &iCols, NULL, NULL, &iCols, NULL, &iCols, &opt, &query,NULL, &info );
+  *optWorkSize= (int)opt.r;
+  *minWorkSize= Max(1, 2*iCols);
+  return info;
+}
+
+static int dsyevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
+{
+  int info=0, query=-1;
+  double opt;
+  C2F(dsyev)("N", "U", &iCols, NULL, &iCols, NULL, &opt, &query, &info);
+  *optWorkSize= (int)opt;
+  *minWorkSize= Max(1, 3*iCols -1);
+  return info;
+}
+
+static int dgeevWorkSizes(int iCols, int lhs, int* optWorkSize, int* minWorkSize)
+{
+  int info=0, query=-1;
+  double opt;
+  C2F(dgeev)("N", "N", &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
+  *optWorkSize= (int)opt;
+  *minWorkSize= (lhs==2) ? Max(1, 4* iCols) : Max(1, 3*iCols);
+  return info;
+}
+
+/**
+ * try to MALLOC optimal (ws[0]) or at least minimal (ws[1]) number of doublecomplex ou double (according to isCplx arg)
+ * @param int const ws[2] in : [0] is the optimal number of elements [1] the minimal
+ * @param int (bool really) isCplx in : if true, allocate for doublecomplex elements, otherwise for double
+ *
+ * @param int* allocated out : nb of allocated elements
+ * @return adress of MALLOCated memory or NULL
+ */
+static void* allocWorkspace(int const ws[2], int const isCplx, int* allocated)
+{
+  int i;
+  void* res=NULL;
+  for(i=0; res==NULL && i!=2; ++i)
+    {
+      res= MALLOC(ws[i] * (isCplx ? sizeof(doublecomplex) : sizeof(double)));
+    }
+  *allocated = (res==NULL ? 0 : ws[i-1]);
+  return res;
+}
+
+/*
+ * internal wrappers around LAPACK function calls
+ * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
+ * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
+ * 
+ * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
+ * @param iCols int nb of rows/cols
+ * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
+
+ * for symetric matrix, eigen values are real
+ * @param pEigenValues double* out ptr to output eigen values
+
+ * for unsymmetric matrix, eigen values are complex :
+ *  in'c' format for unsymmetric real matrix :
+ * @param pEigenValuesReal double* 
+ * @param pEigenValuesImg double*
+ * in 'z' format for unsymmetric complex matrix
+ * @param pEigenValues doublecomplex*
+
+ * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
+
+ * 
+ * @param pWork doublecomplex* scratch ptr to workspace
+ * @param iWorkSize int size of workspace
+ * @param pRWork double* scratch workspace : only used for complex data
+ */
+static int iEigen1CS(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues, doublecomplex* pWork, int iWorkSize, double* pRWork)
+{
+  int info;
+  C2F(zheev)( computeEigenVectors ? "V" : "N", "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, pRWork, &info );
+  return info;
+}
+
+static int iEigen1C(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pRightVectors
+                   , doublecomplex* pWork, int iWorkSize, double* pRWork)
+{
+  int info;
+  C2F(zgeev)( "N", (pRightVectors==NULL ? "N" : "V"), &iCols, pData, &iCols, pEigenValues, NULL, &iCols, pRightVectors, &iCols
+             , pWork, &iWorkSize, pRWork, &info );
+  return info;
+}
+
+static int iEigen1RS(double* pData, int iCols, int computeEigenVectors, double* pEigenValues, double* pWork, int iWorkSize)
+{
+  int info;
+  C2F(dsyev)( (computeEigenVectors ? "V" : "N"), "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, &info );
+  return info;
+}
+
+static int iEigen1R(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pRightVectors
+            , double* pWork, int iWorkSize)
+{
+  int info;
+  C2F(dgeev)( "N", (pRightVectors==NULL ? "N":"V"), &iCols, pData, &iCols, pEigenValuesReal, pEigenValuesImg
+             , NULL, &iCols, pRightVectors, &iCols, pWork, &iWorkSize, &info );
+  return info;
+}
+
+
+/********************************************************************************
+ * functions used by client code to compute eigen values / vectors. See eigen.h *
+ *******************************************************************************/
+
+/*
+ * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complex symmetric matrix.
+ */
+int iEigen1ComplexSymmetricM(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues)
+{
+  int ret=0;
+  int ws[2];
+  int worksize;
+  zheevWorkSizes(iCols, ws, ws+1);
+  {
+    doublecomplex* pWork;
+    double* pRWork;
+    pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
+    pRWork= (double*)MALLOC(Max(1, 3*iCols-2)* sizeof(double));
+    ret= (pWork && pRWork) ? iEigen1CS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize, pRWork ) : 1;
+    FREE(pRWork);
+    FREE(pWork);
+  }
+  return ret;
+}
+/*
+ * Compute eigenvalues (and eigenvectors if pEigenVectors is not NULL) of a complex unsymmetric matrix.
+ */
+int iEigen1ComplexM(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pEigenVectors)
+{
+  int ret=0;
+  int ws[2];
+  int worksize;
+  int lhs= (pEigenVectors == NULL ? 1 :2);
+  zgeevWorkSizes(iCols,lhs, ws, ws+1);
+  {
+    doublecomplex* pWork;
+    double* pRWork;
+    pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
+    pRWork= (double*)MALLOC(2*iCols * sizeof(double));
+    ret = (pWork && pRWork ) ? iEigen1C(pData, iCols, pEigenValues, pEigenVectors, pWork, worksize, pRWork) : 1 ;
+    FREE(pWork);
+    FREE(pRWork);
+  }
+  return ret;
+}
+/*
+ * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complexreal symmetric matrix.
+ */
+int iEigen1RealSymmetricM(double* pData, int iCols, int computeEigenVectors, double* pEigenValues)
+{
+  int ret=0;
+  int ws[2];
+  int worksize;
+  dsyevWorkSizes(iCols, ws, ws+1);
+  {
+    double* pWork;
+    pWork= allocWorkspace(ws, 0, &worksize);
+    ret= (pWork) ? iEigen1RS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize ) : 1 ;
+    FREE(pWork);
+  }
+  return ret;
+}
+/*
+ * Compute eigenvalues (and eigenvectors if pEigenVectorsReal is not NULL) of a real unsymmetric matrix.
+ */
+int iEigen1RealM(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pEigenVectorsReal, double* pEigenVectorsImg)
+{
+  int ret=0;
+  int ws[2];
+  int worksize;
+  int lhs= (pEigenVectorsReal==NULL ? 1 : 2);
+  dgeevWorkSizes(iCols, lhs, ws, ws+1);
+  {
+    double* pWork;
+    double* pRightVectors;
+    pWork= allocWorkspace(ws, 0, &worksize);
+    pRightVectors=  (lhs==2 ? (double*)MALLOC(iCols * iCols * sizeof(double)) : NULL );
+    iEigen1R(pData, iCols, pEigenValuesReal, pEigenValuesImg, pRightVectors, pWork, worksize);
+    if(lhs==2)
+      {
+       assembleEigenvectorsSourceToTarget(iCols, pEigenValuesImg, 
+                                          pRightVectors, 
+                                          pEigenVectorsReal, pEigenVectorsImg);
+       FREE(pRightVectors);
+      }
+  }
+  return ret;
+}
+
+
+
+/*
+ * done in place because a matrix is aways larger than the vector of it's diagonal elements. We must copy from the end to avoid
+ * overwriting elements.
+ * Part of the API: cf eigen.h
+ */
+void expandToDiagonalOfMatrix(double* pData, int iCols)
+{
+  double* ptrDest= pData + iCols*iCols;
+  double const* ptrSrc= pData+iCols;
+  while(--ptrSrc != pData)/* last (first in fact, as we must go backward) diagonal element does not need to be copied as it is already in place*/
+    {
+      *(--ptrDest)= *ptrSrc;
+      ptrDest-= iCols;
+      memset(ptrDest, 0, iCols*sizeof(double));
+    }
+}
+
+/*
+ * complex case is not done inplace because real or imaginary 1x1 matrix are smaller than their complex diagonal.
+ * Part of the API: cf eigen.h
+ */
+void expandZToDiagonalOfCMatrix(doublecomplex const* pZVector, int iCols, double* pRMatrix, double* pIMatrix)
+{
+  double* ptrDestR= pRMatrix + iCols * iCols;
+  double* ptrDestI= pIMatrix + iCols * iCols;
+  double const* ptrSrc= (double const *)(pZVector + iCols);
+  /* we handle the last (in fact first) diagonal element out of the loop because then no bzero is needed */
+  double const* const end= (double const*const)(pZVector+1);
+  while(ptrSrc != end)
+    {
+      *(--ptrDestI)= *(--ptrSrc);
+      *(--ptrDestR)= *(--ptrSrc);
+      ptrDestI-= iCols;
+      bzero(ptrDestI, iCols*sizeof(double));
+      ptrDestR-= iCols;
+      bzero(ptrDestR, iCols*sizeof(double));
+    }
+  /* copy first diagonal element */
+  *pIMatrix= *(--ptrSrc);
+  *pRMatrix= *(--ptrSrc);
+}
+
+/*
+ * Wrappers allocation workspace (after querying optimal worsksize and computing minimal worksize) for LAPACK functions handling two matrix
+ * (contrary to one matrix, there is no sense in sparating query and allocation as it would not improve code factorisation) 
+ *
+ * @param iCols int nb of rows/columns of the matrix
+ *
+ * out:
+ *
+ * @param pWorkSize int* nb of allocated elements
+ *
+ *@return double[complex]* ptr to allocated workspace (NULL if malloc failed)
+ */
+static double* allocateDggevWorkspace(int iCols, int* pWorksize)
+{
+  double* ret=NULL;
+  int info;
+  int query= -1; 
+  double opt;
+  C2F(dggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
+  *pWorksize= (int)opt;
+  ret= MALLOC(*pWorksize);
+  if(!ret)
+    {
+      *pWorksize= Max(1, 8*iCols);
+      ret= MALLOC(*pWorksize);
+      if(!ret)
+       {
+         *pWorksize= 0;
+       }
+    }
+  return ret;
+}
+static doublecomplex* allocateZggevWorkspace(int iCols, int* pWorksize)
+{
+  doublecomplex* ret=NULL;
+  int info;
+  int query= -1; 
+  doublecomplex opt;
+  C2F(zggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, NULL, &info);
+  *pWorksize= (int) opt.r;
+  ret= MALLOC(*pWorksize);
+  if(!ret)
+    {
+      *pWorksize= Max(1, 8*iCols);
+      ret= MALLOC(*pWorksize);
+      if(!ret)
+       {
+         *pWorksize= 0;
+       }
+    }
+  return ret;
+}
+
+/*
+ * internal wrappers around LAPACK function calls
+ * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
+ * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
+ * 
+ * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
+ * @param iCols int nb of rows/cols
+ * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
+
+ * for symetric matrix, eigen values are real
+ * @param pEigenValues double* out ptr to output eigen values
+
+ * for unsymmetric matrix, eigen values are complex :
+ *  in'c' format for unsymmetric real matrix :
+ * @param pEigenValuesReal double* 
+ * @param pEigenValuesImg double*
+ * in 'z' format for unsymmetric complex matrix
+ * @param pEigenValues doublecomplex*
+
+ * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
+
+ * 
+ * @param pWork doublecomplex* scratch ptr to workspace
+ * @param iWorkSize int size of workspace
+ * @param pRWork double* scratch workspace : only used for complex data
+ */
+
+static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
+                  , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork)
+{
+  int info;
+  C2F(zggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlpha, pBeta
+            , pL, &iCols, pR, &iCols, pWork, &iWorkSize, pRWork, &info);
+  return info;
+}
+
+static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
+           , double* pR, double* pL, double* pWork, int iWorkSize)
+{
+  int info;
+  C2F(dggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlphaReal, pAlphaImg
+            , pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, &info);
+  return info;
+}
+
+/******************************************************************************
+ * Part of the API: cf eigen.h
+ ******************************************************************************/
+int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols
+             , doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL)
+{
+  int ret= 0;
+  int iRWorkSize= 0;
+  int worksize;
+  double* pRWork= NULL;
+  doublecomplex* pWork=NULL;  
+  int onlyOneLhs= (pBeta == NULL); /* if beta was not requested (only one lhs), memory was not provided for beta, but will be needed to scale alpha */
+
+  iRWorkSize = Max(1,8*iCols);
+  ret=  ( (pBeta= (onlyOneLhs ? (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex)) : pBeta)) == NULL) 
+    || ( (pRWork = (double*)MALLOC(iRWorkSize * sizeof(double) )) == NULL)
+    || ( (pWork= allocateZggevWorkspace(iCols, &worksize)) == NULL );
+  ret= ret 
+    ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
+    : iEigen2Complex(pData1, pData2, iCols, pAlpha, pBeta, pR, pL, pWork, worksize, pRWork)
+    ;
+  if((ret >=0) && (ret <= iCols) && onlyOneLhs)
+    {/* no error, maybe warnings  and only one lhs : adjust alpha */
+      int ierr; /* original code doe not check ierr */
+      int const delta= sizeof(doublecomplex)/sizeof(double); /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
+      C2F(wwrdiv)(&pAlpha[0].r, &pAlpha[0].i, &delta, &pBeta[0].r, &pBeta[0].i, &delta, &pAlpha[0].r, &pAlpha[0].i, &delta, &iCols, &ierr);
+    }
+  FREE(pRWork);
+  FREE(pWork);
+  if(onlyOneLhs)
+    {
+      FREE(pBeta);
+    }
+  return ret;
+}
+
+/******************************************************************************
+ * Part of the API: cf eigen.h
+ ******************************************************************************/
+int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
+            , double* pRReal, double* pRImg, double* pLReal, double* pLImg)
+{
+  int ret= 0;
+  int worksize;
+  double* pWork=NULL;  
+
+  int onlyOneLhs= (pBeta == NULL);
+  ret=  ( (pBeta= (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
+    || ( (pWork= allocateDggevWorkspace(iCols, &worksize)) == NULL );
+  ret= ret 
+    ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
+    : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize)
+    ;
+  if((ret >=0) && (ret <= iCols) )
+    {/* no error, maybe warnings  and only one lhs : adjust alpha */
+      if(onlyOneLhs)
+       {
+         int i;
+         for(i= 0; i!=iCols; ++i)
+           {
+             pAlphaReal[i] /= pBeta[i];
+             pAlphaImg[i] /= pBeta[i];
+           }
+       }
+      if(pRReal)
+       {
+         assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
+       }
+      if(pLReal)
+       {
+         assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
+       }
+    }
+  FREE(pWork);
+  if(onlyOneLhs)
+    {
+      FREE(pBeta);
+    }
+  return ret;
+}
+
+
+
+/******************************************************************************
+ * Code below lifted from assembleEigenvectors.c 
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2008 - INRIA - MichaĆ«l Baudin
+ *
+ ******************************************************************************/
+//
+// assembleEigenvectorsSourceToTarget --
+//   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
+//   Dispatch the real part of the eigenvectors into the real and imaginary parts,
+//   depending on the imaginary part of the eigenvalues.
+//   The current function reorder the values in the eigenvector array 
+//   and convert from Lapack, compact storage to the Scilab, natural storage.
+// Arguments
+//   iRows : number of rows
+//   eigenvaluesImg, input : the imaginary part of the eigenvalues
+//   EVRealSource, input : the real parts of the source eigenvectors
+//   EVRealTarget, output : real part of the target eigenvectors
+//   EVImgTarget, output : imaginary part of the target eigenvectors
+// Notes
+//   In some eigenvalue computing routines, such as dggev for example, 
+//   the eigenvectors are :
+//   - real (that is with an imaginary part = 0),
+//   - conjugated.
+//   In that case, Lapack stores the eigenvectors in the following compact way.
+//   (Extracted from DGGEV comments)
+//------------------------
+//   The right eigenvectors v(j) are stored one
+//   after another in the columns of VR, in the same order as
+//   their eigenvalues. If the j-th eigenvalue is real, then
+//   v(j) = VR(:,j), the j-th column of VR. If the j-th and
+//   (j+1)-th eigenvalues form a complex conjugate pair, then
+//   v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
+//------------------------
+//   But in Scilab, the eigenvectors must be order in a more natural order,
+//   and this is why a conversion must be performed.
+//
+static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg, 
+                                                                          double * EVRealSource, 
+                                                                          double * EVRealTarget, double * EVImgTarget)
+{      
+
+       double ZERO = 0;
+       int i;
+       int ij;
+       int ij1;
+       int j;
+
+       j = 0;
+       while (j<iRows)
+       {
+               if (eigenvaluesImg[j]==ZERO)
+               {
+                       for(i = 0 ; i < iRows ; i++)
+                       {
+                               ij = i + j * iRows;
+                               EVRealTarget[ij] = EVRealSource[ij];
+                               EVImgTarget[ij] = ZERO;
+                       }
+                       j = j + 1;
+               }
+               else
+               {
+                       for(i = 0 ; i < iRows ; i++)
+                       {
+                               ij = i + j * iRows;
+                               ij1 = i + (j + 1) * iRows;
+                               EVRealTarget[ij] = EVRealSource[ij];
+                               EVImgTarget[ij] = EVRealSource[ij1];
+                               EVRealTarget[ij1] = EVRealSource[ij];
+                               EVImgTarget[ij1] = -EVRealSource[ij1];
+                       }
+                       j = j + 2;
+               }
+       }
+       return 0;
+}
+//
+// assembleEigenvectorsInPlace --
+//   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
+//   Dispatch the real part of the eigenvectors into the real and imaginary parts,
+//   depending on the imaginary part of the eigenvalues.
+//   The current function reorder the values in the eigenvector array 
+//   and convert from Lapack, compact storage to the Scilab, natural storage.
+//   Perform the assembling in place, that is, update the matrix.
+// Arguments
+//   iRows : number of rows
+//   iCols : number of columns
+//   eigenvaluesImg, input : the imaginary part of the eigenvalues
+//   EVReal, input/output : real part of the eigenvectors
+//   EVImg, output : imaginary part of the eigenvectors
+//
+static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
+{      
+
+       double ZERO = 0;
+       int j;
+       int INCY;
+       int totalsize;
+
+       totalsize = iRows * iRows;
+
+       INCY = 1;
+       /*      C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
+       memset(EVImg, 0, totalsize*sizeof(double));
+       j = 0;
+       INCY = 1;
+       while (j<iRows)
+       {
+               if (eigenvaluesImg[j]==ZERO)
+               {
+                       j = j + 1;
+               }
+               else
+               {
+                       int i;
+                       int ij;
+                       int ij1;
+                       for(i = 0 ; i < iRows ; i++)
+                       {
+                               ij = i + j * iRows;
+                               ij1 = i + (j + 1) * iRows;
+                               EVImg[ij]   =   EVReal[ij1];
+                               EVImg[ij1]  = - EVReal[ij1];
+                               EVReal[ij1] =   EVReal[ij];
+                       }
+                       j = j + 2;
+               }
+       }
+       return 0;
+}