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