linear_algebra plugged.
[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  * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE 
5  *
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  */
13 #include <string.h> /* memset */
14 #include "core_math.h"
15 #include "machine.h"
16 #include "MALLOC.h"
17
18 #include "eigen.h"
19
20 /*
21  * Auxiliary functions to expand Lapack compact representation of eigen vectors into scilab matrix
22  * see implementation at the end of the file for API documentation
23  */
24 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg, 
25                                               double * EVRealSource, double * EVRealTarget, double * EVImgTarget);
26 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg);
27
28
29 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
30                    , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork);
31
32 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta, 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
357     C2F(dggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
358
359     *pWorksize= (int)opt;
360     ret= MALLOC(*pWorksize * sizeof(double));
361
362     if(!ret)
363     {
364         *pWorksize = Max(1, 8*iCols);
365         ret= MALLOC(*pWorksize * sizeof(double));
366         if(!ret)
367         {
368             *pWorksize= 0;
369         }
370     }
371     return ret;
372 }
373 static doublecomplex* allocateZggevWorkspace(int iCols, int* pWorksize)
374 {
375     doublecomplex* ret=NULL;
376     int info;
377     int query= -1; 
378     doublecomplex opt;
379
380     C2F(zggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, NULL, &info);
381
382     *pWorksize= (int) opt.r;
383     ret= MALLOC(*pWorksize * sizeof(double));
384
385     if(!ret)
386     {
387         *pWorksize = Max(1, 8*iCols);
388         ret= MALLOC(*pWorksize * sizeof(double));
389         if(!ret)
390         {
391             *pWorksize= 0;
392         }
393     }
394     return ret;
395 }
396
397 /*
398  * internal wrappers around LAPACK function calls
399  * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
400  * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
401  * 
402  * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
403  * @param iCols int nb of rows/cols
404  * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
405
406  * for symetric matrix, eigen values are real
407  * @param pEigenValues double* out ptr to output eigen values
408
409  * for unsymmetric matrix, eigen values are complex :
410  *  in'c' format for unsymmetric real matrix :
411  * @param pEigenValuesReal double* 
412  * @param pEigenValuesImg double*
413  * in 'z' format for unsymmetric complex matrix
414  * @param pEigenValues doublecomplex*
415
416  * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
417
418  * 
419  * @param pWork doublecomplex* scratch ptr to workspace
420  * @param iWorkSize int size of workspace
421  * @param pRWork double* scratch workspace : only used for complex data
422  */
423
424 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork)
425 {
426     int info;
427     C2F(zggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlpha, pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, pRWork, &info);
428     return info;
429 }
430
431 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta, double* pR, double* pL, double* pWork, int iWorkSize)
432 {
433     int info;
434     C2F(dggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlphaReal, pAlphaImg, pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, &info);
435     return info;
436 }
437
438 /******************************************************************************
439  * Part of the API: cf eigen.h
440  ******************************************************************************/
441 int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL)
442 {
443     int ret= 0;
444     int iRWorkSize= 0;
445     int worksize;
446     double* pRWork= NULL;
447     doublecomplex* pWork=NULL;  
448     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 */
449     iRWorkSize = Max(1,8*iCols);
450
451     ret =   ((pBeta = (onlyOneLhs ? (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex)) : pBeta)) == NULL) 
452          || ((pRWork = (double*)MALLOC(iRWorkSize * sizeof(double))) == NULL)
453          || ((pWork = allocateZggevWorkspace(iCols, &worksize)) == NULL);
454
455     /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
456     ret = ret ? -1 : iEigen2Complex(pData1, pData2, iCols, pAlpha, pBeta, pR, pL, pWork, worksize, pRWork);
457
458     if((ret >=0) && (ret <= iCols) && onlyOneLhs)
459     {/* no error, maybe warnings  and only one lhs : adjust alpha */
460         int ierr; /* original code doe not check ierr */
461
462         /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
463         int const delta= sizeof(doublecomplex)/sizeof(double);
464
465         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);
466     }
467
468     FREE(pRWork);
469     FREE(pWork);
470
471     if(onlyOneLhs)
472     {
473         FREE(pBeta);
474     }
475     return ret;
476 }
477
478 /******************************************************************************
479  * Part of the API: cf eigen.h
480  ******************************************************************************/
481 int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta, double* pRReal, double* pRImg, double* pLReal, double* pLImg)
482 {
483     int ret = 0;
484     int worksize;
485     double* pWork = NULL;  
486     int onlyOneLhs = (pBeta == NULL);
487
488     ret =   ((pBeta = (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
489         ||  ((pWork = allocateDggevWorkspace(iCols, &worksize)) == NULL);
490
491     // use explicit -1 to signal malloc error as >0 codes are used for convergence troubles 
492     ret = ret ? -1 : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize);
493
494     if((ret >= 0) && (ret <= iCols))
495     {// no error, maybe warnings  and only one lhs : adjust alpha
496         if(onlyOneLhs)
497         {
498             int i;
499             for(i= 0; i!=iCols; ++i)
500             {
501                 pAlphaReal[i] /= pBeta[i];
502                 pAlphaImg[i] /= pBeta[i];
503             }
504         }
505         if(pRReal)
506         {
507             assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
508         }
509         if(pLReal)
510         {
511             assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
512         }
513     }
514
515     FREE(pWork);
516     if(onlyOneLhs)
517     {
518         FREE(pBeta);
519     }
520     return ret;
521 }
522
523
524
525 /******************************************************************************
526  * Code below lifted from assembleEigenvectors.c 
527  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
528  * Copyright (C) 2008 - INRIA - MichaĆ«l Baudin
529  *
530  ******************************************************************************/
531 //
532 // assembleEigenvectorsSourceToTarget --
533 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
534 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
535 //   depending on the imaginary part of the eigenvalues.
536 //   The current function reorder the values in the eigenvector array 
537 //   and convert from Lapack, compact storage to the Scilab, natural storage.
538 // Arguments
539 //   iRows : number of rows
540 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
541 //   EVRealSource, input : the real parts of the source eigenvectors
542 //   EVRealTarget, output : real part of the target eigenvectors
543 //   EVImgTarget, output : imaginary part of the target eigenvectors
544 // Notes
545 //   In some eigenvalue computing routines, such as dggev for example, 
546 //   the eigenvectors are :
547 //   - real (that is with an imaginary part = 0),
548 //   - conjugated.
549 //   In that case, Lapack stores the eigenvectors in the following compact way.
550 //   (Extracted from DGGEV comments)
551 //------------------------
552 //   The right eigenvectors v(j) are stored one
553 //   after another in the columns of VR, in the same order as
554 //   their eigenvalues. If the j-th eigenvalue is real, then
555 //   v(j) = VR(:,j), the j-th column of VR. If the j-th and
556 //   (j+1)-th eigenvalues form a complex conjugate pair, then
557 //   v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
558 //------------------------
559 //   But in Scilab, the eigenvectors must be order in a more natural order,
560 //   and this is why a conversion must be performed.
561 //
562 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg, 
563                                                                            double * EVRealSource, 
564                                                                            double * EVRealTarget, double * EVImgTarget)
565 {       
566
567         double ZERO = 0;
568         int i;
569         int ij;
570         int ij1;
571         int j;
572
573         j = 0;
574         while (j<iRows)
575         {
576                 if (eigenvaluesImg[j]==ZERO)
577                 {
578                         for(i = 0 ; i < iRows ; i++)
579                         {
580                                 ij = i + j * iRows;
581                                 EVRealTarget[ij] = EVRealSource[ij];
582                                 EVImgTarget[ij] = ZERO;
583                         }
584                         j = j + 1;
585                 }
586                 else
587                 {
588                         for(i = 0 ; i < iRows ; i++)
589                         {
590                                 ij = i + j * iRows;
591                                 ij1 = i + (j + 1) * iRows;
592                                 EVRealTarget[ij] = EVRealSource[ij];
593                                 EVImgTarget[ij] = EVRealSource[ij1];
594                                 EVRealTarget[ij1] = EVRealSource[ij];
595                                 EVImgTarget[ij1] = -EVRealSource[ij1];
596                         }
597                         j = j + 2;
598                 }
599         }
600         return 0;
601 }
602 //
603 // assembleEigenvectorsInPlace --
604 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
605 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
606 //   depending on the imaginary part of the eigenvalues.
607 //   The current function reorder the values in the eigenvector array 
608 //   and convert from Lapack, compact storage to the Scilab, natural storage.
609 //   Perform the assembling in place, that is, update the matrix.
610 // Arguments
611 //   iRows : number of rows
612 //   iCols : number of columns
613 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
614 //   EVReal, input/output : real part of the eigenvectors
615 //   EVImg, output : imaginary part of the eigenvectors
616 //
617 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
618 {       
619
620         double ZERO = 0;
621         int j;
622         int INCY;
623         int totalsize;
624
625         totalsize = iRows * iRows;
626
627         INCY = 1;
628         /*      C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
629         memset(EVImg, 0, totalsize*sizeof(double));
630         j = 0;
631         INCY = 1;
632         while (j<iRows)
633         {
634                 if (eigenvaluesImg[j]==ZERO)
635                 {
636                         j = j + 1;
637                 }
638                 else
639                 {
640                         int i;
641                         int ij;
642                         int ij1;
643                         for(i = 0 ; i < iRows ; i++)
644                         {
645                                 ij = i + j * iRows;
646                                 ij1 = i + (j + 1) * iRows;
647                                 EVImg[ij]   =   EVReal[ij1];
648                                 EVImg[ij1]  = - EVReal[ij1];
649                                 EVReal[ij1] =   EVReal[ij];
650                         }
651                         j = j + 2;
652                 }
653         }
654         return 0;
655 }