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