/*
* 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;
}
/*--------------------------------------------------------------------------*/
--- /dev/null
+/*
+ * 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;
+}
projects / scilab.git / commitdiff