scilab/modules/linear_algebra/src/c/eigen.c

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) 2009 - DIGITEO - Bernard HUGUENEY
   4  *
   5  * This file must be used under the terms of the CeCILL.
   6  * This source file is licensed as described in the file COPYING, which
   7  * you should have received as part of this distribution.  The terms
   8  * are also available at
   9  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  10  *
  11  */
  12 #include <string.h> /* memset */
  13 #include "core_math.h"
  14 #include "machine.h"
  15 #include "MALLOC.h"
  16
  17 #include "eigen.h"
  18
  19 /*
  20  * Auxiliary functions to expand Lapack compact representation of eigen vectors into scilab matrix
  21  * see implementation at the end of the file for API documentation
  22  */
  23 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
  24                                               double * EVRealSource, double * EVRealTarget, double * EVImgTarget);
  25 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg);
  26
  27
  28 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
  29                    , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork);
  30
  31 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
  32                 , double* pR, double* pL, double* pWork, int iWorkSize);
  33
  34 static double* allocateDggevWorkspace(int iCols, int* pWorksize);
  35
  36
  37 /*
  38  * LAPACK Fortran functions used to do the real work for eigen values / vectors of one matrix.
  39  */
  40
  41 extern void C2F(zheev)(char const JOBZ[1]/*'N'|'V'*/, char const UPLO[1]/*'U'|'L'*/, int const* N, doublecomplex* A
  42                        , int const* LDA, double* W, doublecomplex* WORK, int const* LWORK, double* RWORK, int* INFO);
  43
  44 extern void C2F(zgeev)(char const JOBVL[1]/*'N'|'V'*/,char const JOBVR[1]/*'N'|'V'*/, int const* N, doublecomplex* A
  45                        , int const* LDA, doublecomplex* W, doublecomplex* VL, int const* LDVL, doublecomplex* VR, int const* LDVR
  46                        ,doublecomplex* WORK, int const* LWORK, double* RWORK, int* INFO);
  47
  48 extern void C2F(dsyev)( char const JOBZ[1]/*'N'|'V'*/, char const UPLO[1]/*'U'|'L'*/, int const* N, double* A
  49                         , int const* LDA, double* W, double* WORK, int const* LWORK, int* INFO);
  50
  51 extern void C2F(dgeev)(char const JOBVL[1]/*'N'|'V'*/,char const JOBVR[1]/*'N'|'V'*/, int const* N, double* A
  52                        , int const* LDA, double* WR, double* WI, double* VL, int const* LDVL, double* VR, int const* LDVR
  53                        ,double* WORK, int const* LWORK, int* INFO);
  54
  55
  56 /*
  57  * LAPACK Fortran functions used to do the real work for eigen values / vectors of two matrix.
  58  */
  59
  60
  61 extern void C2F(dggev)(char const jobVL[1]/*'N'|'V'*/,char const jobVR[1]/*'N'|'V'*/, int const* n, double* a, int const* ldA
  62                   ,double* b, int const* ldB, double* alphaR, double* alphaI, double* beta, double* vl, int const* ldVl, double* vr, int const* ldVr
  63                   ,double* WORK, int const* LWORK, int* INFO);
  64 extern void C2F(zggev)(char const jobVL[1]/*'N'|'V'*/,char const jobVR[1]/*'N'|'V'*/, int const* n, doublecomplex* a, int const* ldA
  65                   ,doublecomplex* b, int const* ldB, doublecomplex* alpha, doublecomplex* beta, doublecomplex* vl, int const* ldVl
  66                   , doublecomplex* vr, int const* ldVr, doublecomplex* WORK, int const* LWORK, double * RWORK, int* INFO);
  67
  68
  69 /*
  70  * external function to perform division of complex vectors
  71  */
  72 extern int C2F(wwrdiv)(double *ar, double *ai, int const *ia, double *br, double *bi, int const *ib
  73                        , double *rr, double *ri, int const *ir, int const *n, int *ierr);
  74
  75
  76 /*
  77  * Wrappers querying optimal worsksize and computing minimal worksize for LAPACK functions
  78  * @param iCols int nb of rows/columns of the matrix
  79  * @param lhs int nb of expected outputs from the spec scilab function, only when useful to compute optimal worksize
  80  * out:
  81  *
  82  * @param optWorkSize int* points to the optimal worksize
  83  * @param minWorkSize int* points to the minimal worksize
  84  */
  85 static int zheevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
  86 {
  87  int info=0, query=-1;
  88   doublecomplex opt;
  89   C2F(zheev)("N", "U", &iCols, NULL, &iCols, NULL, &opt, &query,NULL, &info );
  90   *optWorkSize= (int) opt.r;
  91   *minWorkSize= Max(1, 2*iCols-1);
  92   return info;
  93 }
  94
  95 static int zgeevWorkSizes(int iCols,  int lhs, int* optWorkSize, int* minWorkSize)
  96 {
  97   int info=0, query=-1;
  98   doublecomplex opt;
  99   C2F(zgeev)("N", (lhs == 1 ? "N" : "V"), &iCols, NULL, &iCols, NULL, NULL, &iCols, NULL, &iCols, &opt, &query,NULL, &info );
 100   *optWorkSize= (int)opt.r;
 101   *minWorkSize= Max(1, 2*iCols);
 102   return info;
 103 }
 104
 105 static int dsyevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
 106 {
 107   int info=0, query=-1;
 108   double opt;
 109   C2F(dsyev)("N", "U", &iCols, NULL, &iCols, NULL, &opt, &query, &info);
 110   *optWorkSize= (int)opt;
 111   *minWorkSize= Max(1, 3*iCols -1);
 112   return info;
 113 }
 114
 115 static int dgeevWorkSizes(int iCols, int lhs, int* optWorkSize, int* minWorkSize)
 116 {
 117   int info=0, query=-1;
 118   double opt;
 119   C2F(dgeev)("N", "N", &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
 120   *optWorkSize= (int)opt;
 121
 122   *minWorkSize= (lhs==2) ? Max(1, 4* iCols) : Max(1, 3*iCols);
 123   return info;
 124 }
 125
 126 /**
 127  * try to MALLOC optimal (ws[0]) or at least minimal (ws[1]) number of doublecomplex ou double (according to isCplx arg)
 128  * @param int const ws[2] in : [0] is the optimal number of elements [1] the minimal
 129  * @param int (bool really) isCplx in : if true, allocate for doublecomplex elements, otherwise for double
 130  *
 131  * @param int* allocated out : nb of allocated elements
 132  * @return adress of MALLOCated memory or NULL
 133  */
 134 static void* allocWorkspace(int const ws[2], int const isCplx, int* allocated)
 135 {
 136   int i;
 137   void* res=NULL;
 138   for(i=0; res==NULL && i!=2; ++i)
 139     {
 140       res= MALLOC(ws[i] * (isCplx ? sizeof(doublecomplex) : sizeof(double)));
 141     }
 142   *allocated = (res==NULL ? 0 : ws[i-1]);
 143   return res;
 144 }
 145
 146 /*
 147  * internal wrappers around LAPACK function calls
 148  * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
 149  * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
 150  *
 151  * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
 152  * @param iCols int nb of rows/cols
 153  * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
 154
 155  * for symetric matrix, eigen values are real
 156  * @param pEigenValues double* out ptr to output eigen values
 157
 158  * for unsymmetric matrix, eigen values are complex :
 159  *  in'c' format for unsymmetric real matrix :
 160  * @param pEigenValuesReal double*
 161  * @param pEigenValuesImg double*
 162  * in 'z' format for unsymmetric complex matrix
 163  * @param pEigenValues doublecomplex*
 164
 165  * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
 166
 167  *
 168  * @param pWork doublecomplex* scratch ptr to workspace
 169  * @param iWorkSize int size of workspace
 170  * @param pRWork double* scratch workspace : only used for complex data
 171  */
 172 static int iEigen1CS(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues, doublecomplex* pWork, int iWorkSize, double* pRWork)
 173 {
 174   int info;
 175   C2F(zheev)( computeEigenVectors ? "V" : "N", "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, pRWork, &info );
 176   return info;
 177 }
 178
 179 static int iEigen1C(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pRightVectors
 180                     , doublecomplex* pWork, int iWorkSize, double* pRWork)
 181 {
 182   int info;
 183   C2F(zgeev)( "N", (pRightVectors==NULL ? "N" : "V"), &iCols, pData, &iCols, pEigenValues, NULL, &iCols, pRightVectors, &iCols
 184               , pWork, &iWorkSize, pRWork, &info );
 185   return info;
 186 }
 187
 188 static int iEigen1RS(double* pData, int iCols, int computeEigenVectors, double* pEigenValues, double* pWork, int iWorkSize)
 189 {
 190   int info;
 191   C2F(dsyev)( (computeEigenVectors ? "V" : "N"), "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, &info );
 192   return info;
 193 }
 194
 195 static int iEigen1R(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pRightVectors
 196              , double* pWork, int iWorkSize)
 197 {
 198   int info;
 199   C2F(dgeev)( "N", (pRightVectors==NULL ? "N":"V"), &iCols, pData, &iCols, pEigenValuesReal, pEigenValuesImg
 200               , NULL, &iCols, pRightVectors, &iCols, pWork, &iWorkSize, &info );
 201   return info;
 202 }
 203
 204
 205 /********************************************************************************
 206  * functions used by client code to compute eigen values / vectors. See eigen.h *
 207  *******************************************************************************/
 208
 209 /*
 210  * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complex symmetric matrix.
 211  */
 212 int iEigen1ComplexSymmetricM(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues)
 213 {
 214   int ret=0;
 215   int ws[2];
 216   int worksize;
 217   zheevWorkSizes(iCols, ws, ws+1);
 218   {
 219     doublecomplex* pWork;
 220     double* pRWork;
 221     pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
 222     pRWork= (double*)MALLOC(Max(1, 3*iCols-2)* sizeof(double));
 223     ret= (pWork && pRWork) ? iEigen1CS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize, pRWork ) : 1;
 224     FREE(pRWork);
 225     FREE(pWork);
 226   }
 227   return ret;
 228 }
 229 /*
 230  * Compute eigenvalues (and eigenvectors if pEigenVectors is not NULL) of a complex unsymmetric matrix.
 231  */
 232 int iEigen1ComplexM(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pEigenVectors)
 233 {
 234   int ret=0;
 235   int ws[2];
 236   int worksize;
 237   int lhs= (pEigenVectors == NULL ? 1 :2);
 238   zgeevWorkSizes(iCols,lhs, ws, ws+1);
 239   {
 240     doublecomplex* pWork;
 241     double* pRWork;
 242     pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
 243     pRWork= (double*)MALLOC(2*iCols * sizeof(double));
 244     ret = (pWork && pRWork ) ? iEigen1C(pData, iCols, pEigenValues, pEigenVectors, pWork, worksize, pRWork) : 1 ;
 245     FREE(pWork);
 246     FREE(pRWork);
 247   }
 248   return ret;
 249 }
 250 /*
 251  * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complexreal symmetric matrix.
 252  */
 253 int iEigen1RealSymmetricM(double* pData, int iCols, int computeEigenVectors, double* pEigenValues)
 254 {
 255   int ret=0;
 256   int ws[2];
 257   int worksize;
 258   dsyevWorkSizes(iCols, ws, ws+1);
 259   {
 260     double* pWork;
 261     pWork= allocWorkspace(ws, 0, &worksize);
 262     ret= (pWork) ? iEigen1RS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize ) : 1 ;
 263     FREE(pWork);
 264   }
 265   return ret;
 266 }
 267 /*
 268  * Compute eigenvalues (and eigenvectors if pEigenVectorsReal is not NULL) of a real unsymmetric matrix.
 269  */
 270 int iEigen1RealM(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pEigenVectorsReal, double* pEigenVectorsImg)
 271 {
 272   int ret=0;
 273   int ws[2];
 274   int worksize;
 275   int lhs= (pEigenVectorsReal==NULL ? 1 : 2);
 276   dgeevWorkSizes(iCols, lhs, ws, ws+1);
 277   {
 278     double* pWork;
 279     double* pRightVectors;
 280     pWork= allocWorkspace(ws, 0, &worksize);
 281     pRightVectors=  (lhs==2 ? (double*)MALLOC(iCols * iCols * sizeof(double)) : NULL );
 282     iEigen1R(pData, iCols, pEigenValuesReal, pEigenValuesImg, pRightVectors, pWork, worksize);
 283     if(lhs==2)
 284       {
 285         assembleEigenvectorsSourceToTarget(iCols, pEigenValuesImg,
 286                                            pRightVectors,
 287                                            pEigenVectorsReal, pEigenVectorsImg);
 288         FREE(pRightVectors);
 289       }
 290   }
 291   return ret;
 292 }
 293
 294
 295
 296 /*
 297  * 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
 298  * overwriting elements.
 299  * Part of the API: cf eigen.h
 300  */
 301 void expandToDiagonalOfMatrix(double* pData, int iCols)
 302 {
 303   double* ptrDest= pData + iCols*iCols;
 304   double const* ptrSrc= pData+iCols;
 305   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*/
 306     {
 307       *(--ptrDest)= *ptrSrc;
 308       ptrDest-= iCols;
 309       memset(ptrDest, 0, iCols*sizeof(double));
 310     }
 311 }
 312
 313 /*
 314  * complex case is not done inplace because real or imaginary 1x1 matrix are smaller than their complex diagonal.
 315  * Part of the API: cf eigen.h
 316  */
 317 void expandZToDiagonalOfCMatrix(doublecomplex const* pZVector, int iCols, double* pRMatrix, double* pIMatrix)
 318 {
 319   double* ptrDestR= pRMatrix + iCols * iCols;
 320   double* ptrDestI= pIMatrix + iCols * iCols;
 321   double const* ptrSrc= (double const *)(pZVector + iCols);
 322   /* we handle the last (in fact first) diagonal element out of the loop because then no memset is needed */
 323   double const* const end= (double const*const)(pZVector+1);
 324   while(ptrSrc != end)
 325     {
 326       *(--ptrDestI)= *(--ptrSrc);
 327       *(--ptrDestR)= *(--ptrSrc);
 328       ptrDestI-= iCols;
 329       memset(ptrDestI, 0, iCols*sizeof(double));
 330       ptrDestR-= iCols;
 331       memset(ptrDestR, 0, iCols*sizeof(double));
 332     }
 333   /* copy first diagonal element */
 334   *pIMatrix= *(--ptrSrc);
 335   *pRMatrix= *(--ptrSrc);
 336 }
 337
 338 /*
 339  * Wrappers allocation workspace (after querying optimal worsksize and computing minimal worksize) for LAPACK functions handling two matrix
 340  * (contrary to one matrix, there is no sense in sparating query and allocation as it would not improve code factorisation)
 341  *
 342  * @param iCols int nb of rows/columns of the matrix
 343  *
 344  * out:
 345  *
 346  * @param pWorkSize int* nb of allocated elements
 347  *
 348  *@return double[complex]* ptr to allocated workspace (NULL if malloc failed)
 349  */
 350 static double* allocateDggevWorkspace(int iCols, int* pWorksize)
 351 {
 352   double* ret=NULL;
 353   int info;
 354   int query= -1;
 355   double opt;
 356   C2F(dggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
 357   *pWorksize= (int)opt;
 358   ret= MALLOC(*pWorksize);
 359   if(!ret)
 360     {
 361       *pWorksize= Max(1, 8*iCols);
 362       ret= MALLOC(*pWorksize);
 363       if(!ret)
 364         {
 365           *pWorksize= 0;
 366         }
 367     }
 368   return ret;
 369 }
 370 static doublecomplex* allocateZggevWorkspace(int iCols, int* pWorksize)
 371 {
 372   doublecomplex* ret=NULL;
 373   int info;
 374   int query= -1;
 375   doublecomplex opt;
 376   C2F(zggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, NULL, &info);
 377   *pWorksize= (int) opt.r;
 378   ret= MALLOC(*pWorksize);
 379   if(!ret)
 380     {
 381       *pWorksize= Max(1, 8*iCols);
 382       ret= MALLOC(*pWorksize);
 383       if(!ret)
 384         {
 385           *pWorksize= 0;
 386         }
 387     }
 388   return ret;
 389 }
 390
 391 /*
 392  * internal wrappers around LAPACK function calls
 393  * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
 394  * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
 395  *
 396  * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
 397  * @param iCols int nb of rows/cols
 398  * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
 399
 400  * for symetric matrix, eigen values are real
 401  * @param pEigenValues double* out ptr to output eigen values
 402
 403  * for unsymmetric matrix, eigen values are complex :
 404  *  in'c' format for unsymmetric real matrix :
 405  * @param pEigenValuesReal double*
 406  * @param pEigenValuesImg double*
 407  * in 'z' format for unsymmetric complex matrix
 408  * @param pEigenValues doublecomplex*
 409
 410  * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
 411
 412  *
 413  * @param pWork doublecomplex* scratch ptr to workspace
 414  * @param iWorkSize int size of workspace
 415  * @param pRWork double* scratch workspace : only used for complex data
 416  */
 417
 418 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
 419                    , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork)
 420 {
 421   int info;
 422   C2F(zggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlpha, pBeta
 423              , pL, &iCols, pR, &iCols, pWork, &iWorkSize, pRWork, &info);
 424   return info;
 425 }
 426
 427 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
 428             , double* pR, double* pL, double* pWork, int iWorkSize)
 429 {
 430   int info;
 431   C2F(dggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlphaReal, pAlphaImg
 432              , pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, &info);
 433   return info;
 434 }
 435
 436 /******************************************************************************
 437  * Part of the API: cf eigen.h
 438  ******************************************************************************/
 439 int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols
 440               , doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL)
 441 {
 442   int ret= 0;
 443   int iRWorkSize= 0;
 444   int worksize;
 445   double* pRWork= NULL;
 446   doublecomplex* pWork=NULL;
 447   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 */
 448
 449   iRWorkSize = Max(1,8*iCols);
 450   ret=  ( (pBeta= (onlyOneLhs ? (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex)) : pBeta)) == NULL)
 451     || ( (pRWork = (double*)MALLOC(iRWorkSize * sizeof(double) )) == NULL)
 452     || ( (pWork= allocateZggevWorkspace(iCols, &worksize)) == NULL );
 453   ret= ret
 454     ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
 455     : iEigen2Complex(pData1, pData2, iCols, pAlpha, pBeta, pR, pL, pWork, worksize, pRWork)
 456     ;
 457   if((ret >=0) && (ret <= iCols) && onlyOneLhs)
 458     {/* no error, maybe warnings  and only one lhs : adjust alpha */
 459       int ierr; /* original code doe not check ierr */
 460       int const delta= sizeof(doublecomplex)/sizeof(double); /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
 461       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);
 462     }
 463   FREE(pRWork);
 464   FREE(pWork);
 465   if(onlyOneLhs)
 466     {
 467       FREE(pBeta);
 468     }
 469   return ret;
 470 }
 471
 472 /******************************************************************************
 473  * Part of the API: cf eigen.h
 474  ******************************************************************************/
 475 int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
 476              , double* pRReal, double* pRImg, double* pLReal, double* pLImg)
 477 {
 478   int ret= 0;
 479   int worksize;
 480   double* pWork=NULL;
 481
 482   int onlyOneLhs= (pBeta == NULL);
 483   ret=  ( (pBeta= (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
 484     || ( (pWork= allocateDggevWorkspace(iCols, &worksize)) == NULL );
 485   ret= ret
 486     ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
 487     : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize)
 488     ;
 489   if((ret >=0) && (ret <= iCols) )
 490     {/* no error, maybe warnings  and only one lhs : adjust alpha */
 491       if(onlyOneLhs)
 492         {
 493           int i;
 494           for(i= 0; i!=iCols; ++i)
 495             {
 496               pAlphaReal[i] /= pBeta[i];
 497               pAlphaImg[i] /= pBeta[i];
 498             }
 499         }
 500       if(pRReal)
 501         {
 502           assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
 503         }
 504       if(pLReal)
 505         {
 506           assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
 507         }
 508     }
 509   FREE(pWork);
 510   if(onlyOneLhs)
 511     {
 512       FREE(pBeta);
 513     }
 514   return ret;
 515 }
 516
 517
 518
 519 /******************************************************************************
 520  * Code below lifted from assembleEigenvectors.c
 521  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 522  * Copyright (C) 2008 - INRIA - Michaël Baudin
 523  *
 524  ******************************************************************************/
 525 //
 526 // assembleEigenvectorsSourceToTarget --
 527 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
 528 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
 529 //   depending on the imaginary part of the eigenvalues.
 530 //   The current function reorder the values in the eigenvector array
 531 //   and convert from Lapack, compact storage to the Scilab, natural storage.
 532 // Arguments
 533 //   iRows : number of rows
 534 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
 535 //   EVRealSource, input : the real parts of the source eigenvectors
 536 //   EVRealTarget, output : real part of the target eigenvectors
 537 //   EVImgTarget, output : imaginary part of the target eigenvectors
 538 // Notes
 539 //   In some eigenvalue computing routines, such as dggev for example,
 540 //   the eigenvectors are :
 541 //   - real (that is with an imaginary part = 0),
 542 //   - conjugated.
 543 //   In that case, Lapack stores the eigenvectors in the following compact way.
 544 //   (Extracted from DGGEV comments)
 545 //------------------------
 546 //   The right eigenvectors v(j) are stored one
 547 //   after another in the columns of VR, in the same order as
 548 //   their eigenvalues. If the j-th eigenvalue is real, then
 549 //   v(j) = VR(:,j), the j-th column of VR. If the j-th and
 550 //   (j+1)-th eigenvalues form a complex conjugate pair, then
 551 //   v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
 552 //------------------------
 553 //   But in Scilab, the eigenvectors must be order in a more natural order,
 554 //   and this is why a conversion must be performed.
 555 //
 556 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
 557                                                                            double * EVRealSource,
 558                                                                            double * EVRealTarget, double * EVImgTarget)
 559 {
 560
 561         double ZERO = 0;
 562         int i;
 563         int ij;
 564         int ij1;
 565         int j;
 566
 567         j = 0;
 568         while (j<iRows)
 569         {
 570                 if (eigenvaluesImg[j]==ZERO)
 571                 {
 572                         for(i = 0 ; i < iRows ; i++)
 573                         {
 574                                 ij = i + j * iRows;
 575                                 EVRealTarget[ij] = EVRealSource[ij];
 576                                 EVImgTarget[ij] = ZERO;
 577                         }
 578                         j = j + 1;
 579                 }
 580                 else
 581                 {
 582                         for(i = 0 ; i < iRows ; i++)
 583                         {
 584                                 ij = i + j * iRows;
 585                                 ij1 = i + (j + 1) * iRows;
 586                                 EVRealTarget[ij] = EVRealSource[ij];
 587                                 EVImgTarget[ij] = EVRealSource[ij1];
 588                                 EVRealTarget[ij1] = EVRealSource[ij];
 589                                 EVImgTarget[ij1] = -EVRealSource[ij1];
 590                         }
 591                         j = j + 2;
 592                 }
 593         }
 594         return 0;
 595 }
 596 //
 597 // assembleEigenvectorsInPlace --
 598 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
 599 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
 600 //   depending on the imaginary part of the eigenvalues.
 601 //   The current function reorder the values in the eigenvector array
 602 //   and convert from Lapack, compact storage to the Scilab, natural storage.
 603 //   Perform the assembling in place, that is, update the matrix.
 604 // Arguments
 605 //   iRows : number of rows
 606 //   iCols : number of columns
 607 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
 608 //   EVReal, input/output : real part of the eigenvectors
 609 //   EVImg, output : imaginary part of the eigenvectors
 610 //
 611 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
 612 {
 613
 614         double ZERO = 0;
 615         int j;
 616         int INCY;
 617         int totalsize;
 618
 619         totalsize = iRows * iRows;
 620
 621         INCY = 1;
 622         /*      C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
 623         memset(EVImg, 0, totalsize*sizeof(double));
 624         j = 0;
 625         INCY = 1;
 626         while (j<iRows)
 627         {
 628                 if (eigenvaluesImg[j]==ZERO)
 629                 {
 630                         j = j + 1;
 631                 }
 632                 else
 633                 {
 634                         int i;
 635                         int ij;
 636                         int ij1;
 637                         for(i = 0 ; i < iRows ; i++)
 638                         {
 639                                 ij = i + j * iRows;
 640                                 ij1 = i + (j + 1) * iRows;
 641                                 EVImg[ij]   =   EVReal[ij1];
 642                                 EVImg[ij1]  = - EVReal[ij1];
 643                                 EVReal[ij1] =   EVReal[ij];
 644                         }
 645                         j = j + 2;
 646                 }
 647         }
 648         return 0;
 649 }