SASification of linear_algebra sci_eigen.c for Scilab spec() function
[scilab.git] / 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 }