2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) ????-2009 - INRIA Bernard Hugueney
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution. The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
12 #include <string.h> /* memset */
13 #include "core_math.h"
20 * Auxiliary functions to expand Lapack compact representation of eigen vectors into scilab matrix
21 * see implementation at the end of the file for API documentation
23 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
24 double * EVRealSource, double * EVRealTarget, double * EVImgTarget);
25 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg);
28 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
29 , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork);
31 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
32 , double* pR, double* pL, double* pWork, int iWorkSize);
34 static double* allocateDggevWorkspace(int iCols, int* pWorksize);
38 * LAPACK Fortran functions used to do the real work for eigen values / vectors of one matrix.
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);
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);
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);
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);
57 * LAPACK Fortran functions used to do the real work for eigen values / vectors of two matrix.
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);
70 * external function to perform division of complex vectors
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);
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
82 * @param optWorkSize int* points to the optimal worksize
83 * @param minWorkSize int* points to the minimal worksize
85 static int zheevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
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);
95 static int zgeevWorkSizes(int iCols, int lhs, int* optWorkSize, int* minWorkSize)
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);
105 static int dsyevWorkSizes(int iCols, int* optWorkSize, int* minWorkSize)
107 int info=0, query=-1;
109 C2F(dsyev)("N", "U", &iCols, NULL, &iCols, NULL, &opt, &query, &info);
110 *optWorkSize= (int)opt;
111 *minWorkSize= Max(1, 3*iCols -1);
115 static int dgeevWorkSizes(int iCols, int lhs, int* optWorkSize, int* minWorkSize)
117 int info=0, query=-1;
119 C2F(dgeev)("N", "N", &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
120 *optWorkSize= (int)opt;
121 *minWorkSize= (lhs==2) ? Max(1, 4* iCols) : Max(1, 3*iCols);
126 * try to MALLOC optimal (ws[0]) or at least minimal (ws[1]) number of doublecomplex ou double (according to isCplx arg)
127 * @param int const ws[2] in : [0] is the optimal number of elements [1] the minimal
128 * @param int (bool really) isCplx in : if true, allocate for doublecomplex elements, otherwise for double
130 * @param int* allocated out : nb of allocated elements
131 * @return adress of MALLOCated memory or NULL
133 static void* allocWorkspace(int const ws[2], int const isCplx, int* allocated)
137 for(i=0; res==NULL && i!=2; ++i)
139 res= MALLOC(ws[i] * (isCplx ? sizeof(doublecomplex) : sizeof(double)));
141 *allocated = (res==NULL ? 0 : ws[i-1]);
146 * internal wrappers around LAPACK function calls
147 * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
148 * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
150 * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
151 * @param iCols int nb of rows/cols
152 * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
154 * for symetric matrix, eigen values are real
155 * @param pEigenValues double* out ptr to output eigen values
157 * for unsymmetric matrix, eigen values are complex :
158 * in'c' format for unsymmetric real matrix :
159 * @param pEigenValuesReal double*
160 * @param pEigenValuesImg double*
161 * in 'z' format for unsymmetric complex matrix
162 * @param pEigenValues doublecomplex*
164 * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
167 * @param pWork doublecomplex* scratch ptr to workspace
168 * @param iWorkSize int size of workspace
169 * @param pRWork double* scratch workspace : only used for complex data
171 static int iEigen1CS(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues, doublecomplex* pWork, int iWorkSize, double* pRWork)
174 C2F(zheev)( computeEigenVectors ? "V" : "N", "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, pRWork, &info );
178 static int iEigen1C(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pRightVectors
179 , doublecomplex* pWork, int iWorkSize, double* pRWork)
182 C2F(zgeev)( "N", (pRightVectors==NULL ? "N" : "V"), &iCols, pData, &iCols, pEigenValues, NULL, &iCols, pRightVectors, &iCols
183 , pWork, &iWorkSize, pRWork, &info );
187 static int iEigen1RS(double* pData, int iCols, int computeEigenVectors, double* pEigenValues, double* pWork, int iWorkSize)
190 C2F(dsyev)( (computeEigenVectors ? "V" : "N"), "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, &info );
194 static int iEigen1R(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pRightVectors
195 , double* pWork, int iWorkSize)
198 C2F(dgeev)( "N", (pRightVectors==NULL ? "N":"V"), &iCols, pData, &iCols, pEigenValuesReal, pEigenValuesImg
199 , NULL, &iCols, pRightVectors, &iCols, pWork, &iWorkSize, &info );
204 /********************************************************************************
205 * functions used by client code to compute eigen values / vectors. See eigen.h *
206 *******************************************************************************/
209 * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complex symmetric matrix.
211 int iEigen1ComplexSymmetricM(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues)
216 zheevWorkSizes(iCols, ws, ws+1);
218 doublecomplex* pWork;
220 pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
221 pRWork= (double*)MALLOC(Max(1, 3*iCols-2)* sizeof(double));
222 ret= (pWork && pRWork) ? iEigen1CS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize, pRWork ) : 1;
229 * Compute eigenvalues (and eigenvectors if pEigenVectors is not NULL) of a complex unsymmetric matrix.
231 int iEigen1ComplexM(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pEigenVectors)
236 int lhs= (pEigenVectors == NULL ? 1 :2);
237 zgeevWorkSizes(iCols,lhs, ws, ws+1);
239 doublecomplex* pWork;
241 pWork= (doublecomplex*)allocWorkspace(ws, 1, &worksize);
242 pRWork= (double*)MALLOC(2*iCols * sizeof(double));
243 ret = (pWork && pRWork ) ? iEigen1C(pData, iCols, pEigenValues, pEigenVectors, pWork, worksize, pRWork) : 1 ;
250 * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complexreal symmetric matrix.
252 int iEigen1RealSymmetricM(double* pData, int iCols, int computeEigenVectors, double* pEigenValues)
257 dsyevWorkSizes(iCols, ws, ws+1);
260 pWork= allocWorkspace(ws, 0, &worksize);
261 ret= (pWork) ? iEigen1RS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize ) : 1 ;
267 * Compute eigenvalues (and eigenvectors if pEigenVectorsReal is not NULL) of a real unsymmetric matrix.
269 int iEigen1RealM(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pEigenVectorsReal, double* pEigenVectorsImg)
274 int lhs= (pEigenVectorsReal==NULL ? 1 : 2);
275 dgeevWorkSizes(iCols, lhs, ws, ws+1);
278 double* pRightVectors;
279 pWork= allocWorkspace(ws, 0, &worksize);
280 pRightVectors= (lhs==2 ? (double*)MALLOC(iCols * iCols * sizeof(double)) : NULL );
281 iEigen1R(pData, iCols, pEigenValuesReal, pEigenValuesImg, pRightVectors, pWork, worksize);
284 assembleEigenvectorsSourceToTarget(iCols, pEigenValuesImg,
286 pEigenVectorsReal, pEigenVectorsImg);
296 * 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
297 * overwriting elements.
298 * Part of the API: cf eigen.h
300 void expandToDiagonalOfMatrix(double* pData, int iCols)
302 double* ptrDest= pData + iCols*iCols;
303 double const* ptrSrc= pData+iCols;
304 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 *(--ptrDest)= *ptrSrc;
308 memset(ptrDest, 0, iCols*sizeof(double));
313 * complex case is not done inplace because real or imaginary 1x1 matrix are smaller than their complex diagonal.
314 * Part of the API: cf eigen.h
316 void expandZToDiagonalOfCMatrix(doublecomplex const* pZVector, int iCols, double* pRMatrix, double* pIMatrix)
318 double* ptrDestR= pRMatrix + iCols * iCols;
319 double* ptrDestI= pIMatrix + iCols * iCols;
320 double const* ptrSrc= (double const *)(pZVector + iCols);
321 /* we handle the last (in fact first) diagonal element out of the loop because then no bzero is needed */
322 double const* const end= (double const*const)(pZVector+1);
325 *(--ptrDestI)= *(--ptrSrc);
326 *(--ptrDestR)= *(--ptrSrc);
328 bzero(ptrDestI, iCols*sizeof(double));
330 bzero(ptrDestR, iCols*sizeof(double));
332 /* copy first diagonal element */
333 *pIMatrix= *(--ptrSrc);
334 *pRMatrix= *(--ptrSrc);
338 * Wrappers allocation workspace (after querying optimal worsksize and computing minimal worksize) for LAPACK functions handling two matrix
339 * (contrary to one matrix, there is no sense in sparating query and allocation as it would not improve code factorisation)
341 * @param iCols int nb of rows/columns of the matrix
345 * @param pWorkSize int* nb of allocated elements
347 *@return double[complex]* ptr to allocated workspace (NULL if malloc failed)
349 static double* allocateDggevWorkspace(int iCols, int* pWorksize)
355 C2F(dggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
356 *pWorksize= (int)opt;
357 ret= MALLOC(*pWorksize);
360 *pWorksize= Max(1, 8*iCols);
361 ret= MALLOC(*pWorksize);
369 static doublecomplex* allocateZggevWorkspace(int iCols, int* pWorksize)
371 doublecomplex* ret=NULL;
375 C2F(zggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, NULL, &info);
376 *pWorksize= (int) opt.r;
377 ret= MALLOC(*pWorksize);
380 *pWorksize= Max(1, 8*iCols);
381 ret= MALLOC(*pWorksize);
391 * internal wrappers around LAPACK function calls
392 * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
393 * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
395 * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
396 * @param iCols int nb of rows/cols
397 * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
399 * for symetric matrix, eigen values are real
400 * @param pEigenValues double* out ptr to output eigen values
402 * for unsymmetric matrix, eigen values are complex :
403 * in'c' format for unsymmetric real matrix :
404 * @param pEigenValuesReal double*
405 * @param pEigenValuesImg double*
406 * in 'z' format for unsymmetric complex matrix
407 * @param pEigenValues doublecomplex*
409 * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
412 * @param pWork doublecomplex* scratch ptr to workspace
413 * @param iWorkSize int size of workspace
414 * @param pRWork double* scratch workspace : only used for complex data
417 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
418 , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork)
421 C2F(zggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlpha, pBeta
422 , pL, &iCols, pR, &iCols, pWork, &iWorkSize, pRWork, &info);
426 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
427 , double* pR, double* pL, double* pWork, int iWorkSize)
430 C2F(dggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlphaReal, pAlphaImg
431 , pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, &info);
435 /******************************************************************************
436 * Part of the API: cf eigen.h
437 ******************************************************************************/
438 int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols
439 , doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL)
444 double* pRWork= NULL;
445 doublecomplex* pWork=NULL;
446 int onlyOneLhs= (pBeta == NULL); /* if beta was not requested (only one lhs), memory was not provided for beta, but will be needed to scale alpha */
448 iRWorkSize = Max(1,8*iCols);
449 ret= ( (pBeta= (onlyOneLhs ? (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex)) : pBeta)) == NULL)
450 || ( (pRWork = (double*)MALLOC(iRWorkSize * sizeof(double) )) == NULL)
451 || ( (pWork= allocateZggevWorkspace(iCols, &worksize)) == NULL );
453 ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
454 : iEigen2Complex(pData1, pData2, iCols, pAlpha, pBeta, pR, pL, pWork, worksize, pRWork)
456 if((ret >=0) && (ret <= iCols) && onlyOneLhs)
457 {/* no error, maybe warnings and only one lhs : adjust alpha */
458 int ierr; /* original code doe not check ierr */
459 int const delta= sizeof(doublecomplex)/sizeof(double); /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
460 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 * Part of the API: cf eigen.h
473 ******************************************************************************/
474 int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
475 , double* pRReal, double* pRImg, double* pLReal, double* pLImg)
481 int onlyOneLhs= (pBeta == NULL);
482 ret= ( (pBeta= (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
483 || ( (pWork= allocateDggevWorkspace(iCols, &worksize)) == NULL );
485 ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
486 : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize)
488 if((ret >=0) && (ret <= iCols) )
489 {/* no error, maybe warnings and only one lhs : adjust alpha */
493 for(i= 0; i!=iCols; ++i)
495 pAlphaReal[i] /= pBeta[i];
496 pAlphaImg[i] /= pBeta[i];
501 assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
505 assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
518 /******************************************************************************
519 * Code below lifted from assembleEigenvectors.c
520 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
521 * Copyright (C) 2008 - INRIA - Michaƫl Baudin
523 ******************************************************************************/
525 // assembleEigenvectorsSourceToTarget --
526 // Assemble conjugated eigenvectors from the real part into real and imaginary parts.
527 // Dispatch the real part of the eigenvectors into the real and imaginary parts,
528 // depending on the imaginary part of the eigenvalues.
529 // The current function reorder the values in the eigenvector array
530 // and convert from Lapack, compact storage to the Scilab, natural storage.
532 // iRows : number of rows
533 // eigenvaluesImg, input : the imaginary part of the eigenvalues
534 // EVRealSource, input : the real parts of the source eigenvectors
535 // EVRealTarget, output : real part of the target eigenvectors
536 // EVImgTarget, output : imaginary part of the target eigenvectors
538 // In some eigenvalue computing routines, such as dggev for example,
539 // the eigenvectors are :
540 // - real (that is with an imaginary part = 0),
542 // In that case, Lapack stores the eigenvectors in the following compact way.
543 // (Extracted from DGGEV comments)
544 //------------------------
545 // The right eigenvectors v(j) are stored one
546 // after another in the columns of VR, in the same order as
547 // their eigenvalues. If the j-th eigenvalue is real, then
548 // v(j) = VR(:,j), the j-th column of VR. If the j-th and
549 // (j+1)-th eigenvalues form a complex conjugate pair, then
550 // v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
551 //------------------------
552 // But in Scilab, the eigenvectors must be order in a more natural order,
553 // and this is why a conversion must be performed.
555 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
556 double * EVRealSource,
557 double * EVRealTarget, double * EVImgTarget)
569 if (eigenvaluesImg[j]==ZERO)
571 for(i = 0 ; i < iRows ; i++)
574 EVRealTarget[ij] = EVRealSource[ij];
575 EVImgTarget[ij] = ZERO;
581 for(i = 0 ; i < iRows ; i++)
584 ij1 = i + (j + 1) * iRows;
585 EVRealTarget[ij] = EVRealSource[ij];
586 EVImgTarget[ij] = EVRealSource[ij1];
587 EVRealTarget[ij1] = EVRealSource[ij];
588 EVImgTarget[ij1] = -EVRealSource[ij1];
596 // assembleEigenvectorsInPlace --
597 // Assemble conjugated eigenvectors from the real part into real and imaginary parts.
598 // Dispatch the real part of the eigenvectors into the real and imaginary parts,
599 // depending on the imaginary part of the eigenvalues.
600 // The current function reorder the values in the eigenvector array
601 // and convert from Lapack, compact storage to the Scilab, natural storage.
602 // Perform the assembling in place, that is, update the matrix.
604 // iRows : number of rows
605 // iCols : number of columns
606 // eigenvaluesImg, input : the imaginary part of the eigenvalues
607 // EVReal, input/output : real part of the eigenvectors
608 // EVImg, output : imaginary part of the eigenvectors
610 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
618 totalsize = iRows * iRows;
621 /* C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
622 memset(EVImg, 0, totalsize*sizeof(double));
627 if (eigenvaluesImg[j]==ZERO)
636 for(i = 0 ; i < iRows ; i++)
639 ij1 = i + (j + 1) * iRows;
640 EVImg[ij] = EVReal[ij1];
641 EVImg[ij1] = - EVReal[ij1];
642 EVReal[ij1] = EVReal[ij];