fix copyrights string format
[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 - 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 }