2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2009 - DIGITEO - 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;
122 *minWorkSize= (lhs==2) ? Max(1, 4* iCols) : Max(1, 3*iCols);
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
131 * @param int* allocated out : nb of allocated elements
132 * @return adress of MALLOCated memory or NULL
134 static void* allocWorkspace(int const ws[2], int const isCplx, int* allocated)
138 for(i=0; res==NULL && i!=2; ++i)
140 res= MALLOC(ws[i] * (isCplx ? sizeof(doublecomplex) : sizeof(double)));
142 *allocated = (res==NULL ? 0 : ws[i-1]);
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
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)
155 * for symetric matrix, eigen values are real
156 * @param pEigenValues double* out ptr to output eigen values
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*
165 * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
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
172 static int iEigen1CS(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues, doublecomplex* pWork, int iWorkSize, double* pRWork)
175 C2F(zheev)( computeEigenVectors ? "V" : "N", "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, pRWork, &info );
179 static int iEigen1C(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pRightVectors
180 , doublecomplex* pWork, int iWorkSize, double* pRWork)
183 C2F(zgeev)( "N", (pRightVectors==NULL ? "N" : "V"), &iCols, pData, &iCols, pEigenValues, NULL, &iCols, pRightVectors, &iCols
184 , pWork, &iWorkSize, pRWork, &info );
188 static int iEigen1RS(double* pData, int iCols, int computeEigenVectors, double* pEigenValues, double* pWork, int iWorkSize)
191 C2F(dsyev)( (computeEigenVectors ? "V" : "N"), "U", &iCols, pData, &iCols, pEigenValues, pWork, &iWorkSize, &info );
195 static int iEigen1R(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pRightVectors
196 , double* pWork, int iWorkSize)
199 C2F(dgeev)( "N", (pRightVectors==NULL ? "N":"V"), &iCols, pData, &iCols, pEigenValuesReal, pEigenValuesImg
200 , NULL, &iCols, pRightVectors, &iCols, pWork, &iWorkSize, &info );
205 /********************************************************************************
206 * functions used by client code to compute eigen values / vectors. See eigen.h *
207 *******************************************************************************/
210 * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complex symmetric matrix.
212 int iEigen1ComplexSymmetricM(doublecomplex* pData, int iCols, int computeEigenVectors, double* pEigenValues)
217 zheevWorkSizes(iCols, ws, ws+1);
219 doublecomplex* pWork;
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;
230 * Compute eigenvalues (and eigenvectors if pEigenVectors is not NULL) of a complex unsymmetric matrix.
232 int iEigen1ComplexM(doublecomplex* pData, int iCols, doublecomplex* pEigenValues, doublecomplex* pEigenVectors)
237 int lhs= (pEigenVectors == NULL ? 1 :2);
238 zgeevWorkSizes(iCols,lhs, ws, ws+1);
240 doublecomplex* pWork;
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 ;
251 * Compute eigenvalues (and eigenvectors if computeEigenVectors is true) of a complexreal symmetric matrix.
253 int iEigen1RealSymmetricM(double* pData, int iCols, int computeEigenVectors, double* pEigenValues)
258 dsyevWorkSizes(iCols, ws, ws+1);
261 pWork= allocWorkspace(ws, 0, &worksize);
262 ret= (pWork) ? iEigen1RS(pData, iCols, computeEigenVectors, pEigenValues, pWork, worksize ) : 1 ;
268 * Compute eigenvalues (and eigenvectors if pEigenVectorsReal is not NULL) of a real unsymmetric matrix.
270 int iEigen1RealM(double* pData, int iCols, double* pEigenValuesReal, double* pEigenValuesImg, double* pEigenVectorsReal, double* pEigenVectorsImg)
275 int lhs= (pEigenVectorsReal==NULL ? 1 : 2);
276 dgeevWorkSizes(iCols, lhs, ws, ws+1);
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);
285 assembleEigenvectorsSourceToTarget(iCols, pEigenValuesImg,
287 pEigenVectorsReal, pEigenVectorsImg);
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
301 void expandToDiagonalOfMatrix(double* pData, int iCols)
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*/
307 *(--ptrDest)= *ptrSrc;
309 memset(ptrDest, 0, iCols*sizeof(double));
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
317 void expandZToDiagonalOfCMatrix(doublecomplex const* pZVector, int iCols, double* pRMatrix, double* pIMatrix)
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);
326 *(--ptrDestI)= *(--ptrSrc);
327 *(--ptrDestR)= *(--ptrSrc);
329 memset(ptrDestI, 0, iCols*sizeof(double));
331 memset(ptrDestR, 0, iCols*sizeof(double));
333 /* copy first diagonal element */
334 *pIMatrix= *(--ptrSrc);
335 *pRMatrix= *(--ptrSrc);
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)
342 * @param iCols int nb of rows/columns of the matrix
346 * @param pWorkSize int* nb of allocated elements
348 *@return double[complex]* ptr to allocated workspace (NULL if malloc failed)
350 static double* allocateDggevWorkspace(int iCols, int* pWorksize)
356 C2F(dggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, &info);
357 *pWorksize= (int)opt;
358 ret= MALLOC(*pWorksize);
361 *pWorksize= Max(1, 8*iCols);
362 ret= MALLOC(*pWorksize);
370 static doublecomplex* allocateZggevWorkspace(int iCols, int* pWorksize)
372 doublecomplex* ret=NULL;
376 C2F(zggev)("N", "N", &iCols, NULL, &iCols, NULL, &iCols, NULL, NULL, NULL, &iCols, NULL, &iCols, &opt, &query, NULL, &info);
377 *pWorksize= (int) opt.r;
378 ret= MALLOC(*pWorksize);
381 *pWorksize= Max(1, 8*iCols);
382 ret= MALLOC(*pWorksize);
392 * internal wrappers around LAPACK function calls
393 * For symmetric cases, we use use an int (bool really) computeEigenVectors to express wether eigenvalues should be computed.
394 * For unsymmetric cases, eigenvectors are not computed in place, so a NULL pointer for ouput indicates that eigen vectors should not be computed
396 * @param pData double[complex]* ptr to data matrix for symmetric matrix, also used as output for eigen vectors (when computed)
397 * @param iCols int nb of rows/cols
398 * @param computeEigenVectors int (boolean semantics) weither eigen vectors should be computed only used for symmetric data (cf. supra)
400 * for symetric matrix, eigen values are real
401 * @param pEigenValues double* out ptr to output eigen values
403 * for unsymmetric matrix, eigen values are complex :
404 * in'c' format for unsymmetric real matrix :
405 * @param pEigenValuesReal double*
406 * @param pEigenValuesImg double*
407 * in 'z' format for unsymmetric complex matrix
408 * @param pEigenValues doublecomplex*
410 * @param pRightVectors double[complex]* output eigen vectors (for unsymmetric matrix), when NULL, eigen vectors are not computed
413 * @param pWork doublecomplex* scratch ptr to workspace
414 * @param iWorkSize int size of workspace
415 * @param pRWork double* scratch workspace : only used for complex data
418 static int iEigen2Complex(doublecomplex* pData1, doublecomplex* pData2, int iCols, doublecomplex* pAlpha, doublecomplex* pBeta
419 , doublecomplex* pR, doublecomplex* pL, doublecomplex* pWork, int iWorkSize, double* pRWork)
422 C2F(zggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlpha, pBeta
423 , pL, &iCols, pR, &iCols, pWork, &iWorkSize, pRWork, &info);
427 static int iEigen2Real(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
428 , double* pR, double* pL, double* pWork, int iWorkSize)
431 C2F(dggev)((pL != NULL) ? "V" : "N", (pR != NULL) ? "V" : "N", &iCols, pData1, &iCols, pData2, &iCols, pAlphaReal, pAlphaImg
432 , pBeta, pL, &iCols, pR, &iCols, pWork, &iWorkSize, &info);
436 /******************************************************************************
437 * Part of the API: cf eigen.h
438 ******************************************************************************/
439 int iEigen2ComplexM(doublecomplex* pData1, doublecomplex* pData2, int iCols
440 , doublecomplex* pAlpha, doublecomplex* pBeta, doublecomplex* pR, doublecomplex* pL)
445 double* pRWork= NULL;
446 doublecomplex* pWork=NULL;
447 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 ret= ( (pBeta= (onlyOneLhs ? (doublecomplex*)MALLOC(iCols * sizeof(doublecomplex)) : pBeta)) == NULL)
451 || ( (pRWork = (double*)MALLOC(iRWorkSize * sizeof(double) )) == NULL)
452 || ( (pWork= allocateZggevWorkspace(iCols, &worksize)) == NULL );
454 ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
455 : iEigen2Complex(pData1, pData2, iCols, pAlpha, pBeta, pR, pL, pWork, worksize, pRWork)
457 if((ret >=0) && (ret <= iCols) && onlyOneLhs)
458 {/* no error, maybe warnings and only one lhs : adjust alpha */
459 int ierr; /* original code doe not check ierr */
460 int const delta= sizeof(doublecomplex)/sizeof(double); /* in fact it's (&pAlpha[1].r - &pAlpha[1].r) / sizeof(double), or... 2 :) */
461 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 * Part of the API: cf eigen.h
474 ******************************************************************************/
475 int iEigen2RealM(double* pData1, double* pData2, int iCols, double* pAlphaReal, double* pAlphaImg, double* pBeta
476 , double* pRReal, double* pRImg, double* pLReal, double* pLImg)
482 int onlyOneLhs= (pBeta == NULL);
483 ret= ( (pBeta= (onlyOneLhs ? (double*)MALLOC(iCols * sizeof(double)) : pBeta)) == NULL)
484 || ( (pWork= allocateDggevWorkspace(iCols, &worksize)) == NULL );
486 ? -1 /* use explicit -1 to signal malloc error as >0 codes are used for convergence troubles */
487 : iEigen2Real(pData1, pData2, iCols, pAlphaReal, pAlphaImg, pBeta, pRReal, pLReal, pWork, worksize)
489 if((ret >=0) && (ret <= iCols) )
490 {/* no error, maybe warnings and only one lhs : adjust alpha */
494 for(i= 0; i!=iCols; ++i)
496 pAlphaReal[i] /= pBeta[i];
497 pAlphaImg[i] /= pBeta[i];
502 assembleEigenvectorsInPlace(iCols, pAlphaImg, pRReal, pRImg);
506 assembleEigenvectorsInPlace(iCols, pAlphaImg, pLReal, pLImg);
519 /******************************************************************************
520 * Code below lifted from assembleEigenvectors.c
521 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
522 * Copyright (C) 2008 - INRIA - Michaƫl Baudin
524 ******************************************************************************/
526 // assembleEigenvectorsSourceToTarget --
527 // Assemble conjugated eigenvectors from the real part into real and imaginary parts.
528 // Dispatch the real part of the eigenvectors into the real and imaginary parts,
529 // depending on the imaginary part of the eigenvalues.
530 // The current function reorder the values in the eigenvector array
531 // and convert from Lapack, compact storage to the Scilab, natural storage.
533 // iRows : number of rows
534 // eigenvaluesImg, input : the imaginary part of the eigenvalues
535 // EVRealSource, input : the real parts of the source eigenvectors
536 // EVRealTarget, output : real part of the target eigenvectors
537 // EVImgTarget, output : imaginary part of the target eigenvectors
539 // In some eigenvalue computing routines, such as dggev for example,
540 // the eigenvectors are :
541 // - real (that is with an imaginary part = 0),
543 // In that case, Lapack stores the eigenvectors in the following compact way.
544 // (Extracted from DGGEV comments)
545 //------------------------
546 // The right eigenvectors v(j) are stored one
547 // after another in the columns of VR, in the same order as
548 // their eigenvalues. If the j-th eigenvalue is real, then
549 // v(j) = VR(:,j), the j-th column of VR. If the j-th and
550 // (j+1)-th eigenvalues form a complex conjugate pair, then
551 // v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
552 //------------------------
553 // But in Scilab, the eigenvectors must be order in a more natural order,
554 // and this is why a conversion must be performed.
556 static int assembleEigenvectorsSourceToTarget(int iRows, double * eigenvaluesImg,
557 double * EVRealSource,
558 double * EVRealTarget, double * EVImgTarget)
570 if (eigenvaluesImg[j]==ZERO)
572 for(i = 0 ; i < iRows ; i++)
575 EVRealTarget[ij] = EVRealSource[ij];
576 EVImgTarget[ij] = ZERO;
582 for(i = 0 ; i < iRows ; i++)
585 ij1 = i + (j + 1) * iRows;
586 EVRealTarget[ij] = EVRealSource[ij];
587 EVImgTarget[ij] = EVRealSource[ij1];
588 EVRealTarget[ij1] = EVRealSource[ij];
589 EVImgTarget[ij1] = -EVRealSource[ij1];
597 // assembleEigenvectorsInPlace --
598 // Assemble conjugated eigenvectors from the real part into real and imaginary parts.
599 // Dispatch the real part of the eigenvectors into the real and imaginary parts,
600 // depending on the imaginary part of the eigenvalues.
601 // The current function reorder the values in the eigenvector array
602 // and convert from Lapack, compact storage to the Scilab, natural storage.
603 // Perform the assembling in place, that is, update the matrix.
605 // iRows : number of rows
606 // iCols : number of columns
607 // eigenvaluesImg, input : the imaginary part of the eigenvalues
608 // EVReal, input/output : real part of the eigenvectors
609 // EVImg, output : imaginary part of the eigenvectors
611 static int assembleEigenvectorsInPlace(int iRows, double * eigenvaluesImg, double * EVReal, double * EVImg)
619 totalsize = iRows * iRows;
622 /* C2F(dset)(&totalsize,&ZERO,EVImg,&INCY);*/
623 memset(EVImg, 0, totalsize*sizeof(double));
628 if (eigenvaluesImg[j]==ZERO)
637 for(i = 0 ; i < iRows ; i++)
640 ij1 = i + (j + 1) * iRows;
641 EVImg[ij] = EVReal[ij1];
642 EVImg[ij1] = - EVReal[ij1];
643 EVReal[ij1] = EVReal[ij];