Massive indent of all codes:
[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     {
460         /* no error, maybe warnings  and only one lhs : adjust alpha */
461         int ierr; /* original code doe not check ierr */
462
463         /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
464         int const delta = sizeof(doublecomplex) / sizeof(double);
465
466         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);
467     }
468
469     FREE(pRWork);
470     FREE(pWork);
471
472     if (onlyOneLhs)
473     {
474         FREE(pBeta);
475     }
476     return ret;
477 }
478
479 /******************************************************************************
480  * Part of the API: cf eigen.h
481  ******************************************************************************/
482 int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta, double* pRReal, double* pRImg, double* pLReal, double* pLImg)
483 {
484     int ret = 0;
485     int worksize;
486     double* pWork = NULL;
487     int onlyOneLhs = (pBeta == NULL);
488
489     ret =   ((pBeta = (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
490             ||  ((pWork = allocateDggevWorkspace(iCols, &worksize)) == NULL);
491
492     // use explicit -1 to signal malloc error as >0 codes are used for convergence troubles
493     ret = ret ? -1 : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize);
494
495     if ((ret >= 0) && (ret <= iCols))
496     {
497         // no error, maybe warnings  and only one lhs : adjust alpha
498         if (onlyOneLhs)
499         {
500             int i;
501             for (i = 0; i != iCols; ++i)
502             {
503                 pAlphaReal[i] /= pBeta[i];
504                 pAlphaImg[i] /= pBeta[i];
505             }
506         }
507         if (pRReal)
508         {
509             assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
510         }
511         if (pLReal)
512         {
513             assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
514         }
515     }
516
517     FREE(pWork);
518     if (onlyOneLhs)
519     {
520         FREE(pBeta);
521     }
522     return ret;
523 }
524
525
526
527 /******************************************************************************
528  * Code below lifted from assembleEigenvectors.c
529  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
530  * Copyright (C) 2008 - INRIA - MichaĆ«l Baudin
531  *
532  ******************************************************************************/
533 //
534 // assembleEigenvectorsSourceToTarget --
535 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
536 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
537 //   depending on the imaginary part of the eigenvalues.
538 //   The current function reorder the values in the eigenvector array
539 //   and convert from Lapack, compact storage to the Scilab, natural storage.
540 // Arguments
541 //   iRows : number of rows
542 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
543 //   EVRealSource, input : the real parts of the source eigenvectors
544 //   EVRealTarget, output : real part of the target eigenvectors
545 //   EVImgTarget, output : imaginary part of the target eigenvectors
546 // Notes
547 //   In some eigenvalue computing routines, such as dggev for example,
548 //   the eigenvectors are :
549 //   - real (that is with an imaginary part = 0),
550 //   - conjugated.
551 //   In that case, Lapack stores the eigenvectors in the following compact way.
552 //   (Extracted from DGGEV comments)
553 //------------------------
554 //   The right eigenvectors v(j) are stored one
555 //   after another in the columns of VR, in the same order as
556 //   their eigenvalues. If the j-th eigenvalue is real, then
557 //   v(j) = VR(:,j), the j-th column of VR. If the j-th and
558 //   (j+1)-th eigenvalues form a complex conjugate pair, then
559 //   v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
560 //------------------------
561 //   But in Scilab, the eigenvectors must be order in a more natural order,
562 //   and this is why a conversion must be performed.
563 //
564 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
565         double * EVRealSource,
566         double * EVRealTarget, double * EVImgTarget)
567 {
568
569     double ZERO = 0;
570     int i;
571     int ij;
572     int ij1;
573     int j;
574
575     j = 0;
576     while (j < iRows)
577     {
578         if (eigenvaluesImg[j] == ZERO)
579         {
580             for (i = 0 ; i < iRows ; i++)
581             {
582                 ij = i + j * iRows;
583                 EVRealTarget[ij] = EVRealSource[ij];
584                 EVImgTarget[ij] = ZERO;
585             }
586             j = j + 1;
587         }
588         else
589         {
590             for (i = 0 ; i < iRows ; i++)
591             {
592                 ij = i + j * iRows;
593                 ij1 = i + (j + 1) * iRows;
594                 EVRealTarget[ij] = EVRealSource[ij];
595                 EVImgTarget[ij] = EVRealSource[ij1];
596                 EVRealTarget[ij1] = EVRealSource[ij];
597                 EVImgTarget[ij1] = -EVRealSource[ij1];
598             }
599             j = j + 2;
600         }
601     }
602     return 0;
603 }
604 //
605 // assembleEigenvectorsInPlace --
606 //   Assemble conjugated eigenvectors from the real part into real and imaginary parts.
607 //   Dispatch the real part of the eigenvectors into the real and imaginary parts,
608 //   depending on the imaginary part of the eigenvalues.
609 //   The current function reorder the values in the eigenvector array
610 //   and convert from Lapack, compact storage to the Scilab, natural storage.
611 //   Perform the assembling in place, that is, update the matrix.
612 // Arguments
613 //   iRows : number of rows
614 //   iCols : number of columns
615 //   eigenvaluesImg, input : the imaginary part of the eigenvalues
616 //   EVReal, input/output : real part of the eigenvectors
617 //   EVImg, output : imaginary part of the eigenvectors
618 //
619 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
620 {
621
622     double ZERO = 0;
623     int j;
624     int INCY;
625     int totalsize;
626
627     totalsize = iRows * iRows;
628
629     INCY = 1;
630     /*  C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
631     memset(EVImg, 0, totalsize * sizeof(double));
632     j = 0;
633     INCY = 1;
634     while (j < iRows)
635     {
636         if (eigenvaluesImg[j] == ZERO)
637         {
638             j = j + 1;
639         }
640         else
641         {
642             int i;
643             int ij;
644             int ij1;
645             for (i = 0 ; i < iRows ; i++)
646             {
647                 ij = i + j * iRows;
648                 ij1 = i + (j + 1) * iRows;
649                 EVImg[ij]   =   EVReal[ij1];
650                 EVImg[ij1]  = - EVReal[ij1];
651                 EVReal[ij1] =   EVReal[ij];
652             }
653             j = j + 2;
654         }
655     }
656     return 0;
657 }