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