src/Matrix/h_matrix.cc
/*************************************************************************
*									**
** The routine diag combines the routines cred,	rred and tqli (using	**
** the function sign) to diagonalize a complex Hermitian matrix.  	**
**									**
** During the course of the transformations, exclusively unitary	**
** (orthogonal) transformations are used. Errors come mainly from	**
** rounding. The resulting eigenbase of a Hermitian array deviates 	**
** from unitarity only negligibly.					**
**									**
** This is THE Entry Point to Hermitian Array Diagonalizations!		**
** This Is The Workhorse Behind Many GAMMA Calculations............!	**
** This Routine Needs to Be Working Very Well......................!	**
**									**
** Note That There Is Currently One Unique Feature Of This Function	**
** That Sets It Apart From Most Other GAMMA Matrix Functions. It Sets	**
** The Pointers mxd and mxev To Point To A Diagonal And A Normal Array	**
** Respectively.  In Turn, The Function Calls Diagonalization Routines	**
** That EXCLUSIVELY TAKE d_matrix and n_matrix!				**
**									**
*************************************************************************/
The routine diag combines the routines cred,   rred and tqli (using
the function sign) to diagonalize a complex Hermitian matrix.

During the course of the transformations, exclusively unitary
(orthogonal) transformations are used. Errors come mainly from
rounding. The resulting eigenbase of a Hermitian array deviates
from unitarity only negligibly.

This is THE Entry Point to Hermitian Array Diagonalizations.
This Is The Workhorse Behind Many GAMMA Calculations.............

Note That There Is Currently One Unique Feature Of This Function
That Sets It Apart From Most Other GAMMA Matrix Functions. It Sets
The Pointers mxd and mxev To Point To A Diagonal And A Normal Array
Respectively.  In Turn, The Function Calls Diagonalization Routines
That EXCLUSIVELY TAKE d_matrix and n_matrix!
The routine diag combines the routines cred rred and tqli ( using the function sign ) to diagonalize a complex Hermitian matrix . During the course of the transformations exclusively unitary ( orthogonal ) transformations are used . Errors come mainly from rounding . The resulting eigenbase of a Hermitian array deviates from unitarity only negligibly . This is THE Entry Point to Hermitian Array Diagonalizations . This Is The Workhorse Behind Many GAMMA Calculations ............. Note That There Is Currently One Unique Feature Of This Function That Sets It Apart From Most Other GAMMA Matrix Functions . It Sets The Pointers mxd and mxev To Point To A Diagonal And A Normal Array Respectively . In Turn The Function Calls Diagonalization Routines That EXCLUSIVELY TAKE d_matrix and n_matrix!
void h_matrix::diag(_matrix* (&D), _matrix* (&U))
 Output       void    : The matrices D and U are set to the diagonal
                          matrix of hmx eigenvalues and the matrix
                          of dmx eigenvectors respectively.
 Note 1               : Both U and D MUST be input pointing to
                          NULL matrix. Their data will be allocated in
                          cred and tqli respectively. Their reference
                          counting MUST be handled externally.
 Note 2               : D will be a d_matrix, U n_matrix unitary
{
#if defined(_USING_SUNPERFLIB_) || defined(_USING_LAPACK_)
    int nrows = this->rows();
    int ncols = this->cols();
    if(nrows < 32)   then do it the old fashioned gamma way.
    {
        h_matrix H(*this);			     Make a copy of hmx to use
        H.cred(U);				         To Hermitian tridiagonal form
        H.rred(U);				         To real symm. tridiagonal form
        H.tqli(U, D);				     To real diagonal form
        ((n_matrix*)U)->unitary = true;	 This should be unitary now
	std::cerr << "h_matrix diag: using gamma code\n";
    }
    else
    {
	n_matrix* hmx =	new n_matrix(nrows, ncols);
	this->convert(hmx);
  	n_matrix* hmx1 = (n_matrix *)hmx->transpose();
#ifdef _USING_SUNPERFLIB_
        char jobz = 'V';   Calculate "eigenvectors" AND "eigenvalues".
        char uplo = 'U';   Upper triangular...
        int  N    = nrows;
        int  lda  = nrows;
        double *w_eig = new double[N];
        int info = -55555;
        zheev(jobz, uplo, N, (doublecomplex *)hmx1->data, lda, w_eig, &info);
#endif
#ifdef _USING_LAPACK_
        char jobz = 'V';   Calculate "eigenvectors" AND "eigenvalues".
        char uplo = 'U';   Upper triangular...
#if defined(__LP64__) 
        int  N    = nrows;
        int  lda  = nrows;
        int info = -55555;
        int lwork = 2*N + 10;
#else
        long int  N    = nrows;
        long int  lda  = nrows;
        long int info = -55555;
        long int lwork = 2*N + 10;
#endif
        double *w_eig = new double[N];
        complex *work= new complex [2*lwork];
        double *rwork = new double [3*N-2];
        zheev_(&jobz, &uplo, &N, (__CLPK_doublecomplex *) hmx1->data, &lda, w_eig, (__CLPK_doublecomplex *) work, &lwork, rwork, &info);
#endif
        if(info == 0)
        {  all is well.
        }
        else if (info > 0)
        { std::cerr << "\nDiagonalization failed to converge\n";
        }
        else if(info == -55555)
        { std::cerr << "\nReturn value, 'info', does not appear to have been set\n";
        }
        else
        { std::cerr << "\ninfo = " << info << "\n";
        }
         Create the matrixes D and U.
        D = new d_matrix(nrows, ncols);
        U = new n_matrix(nrows, ncols);
         Copy results back into D and U.
	complex *hmxp = hmx1->data;
        for(int i=0; i<nrows; i++)
        { D->put( w_eig[i], i, i);
          for(int j=0; j<ncols; j++)
          { U->put(hmxp[j*nrows+i], i, j);
          }
        }
        ((n_matrix*)U)->unitary = true;	 This should be unitary now
	delete hmx;
	delete hmx1;
	delete [] w_eig;
#ifdef _USING_LAPACK_
	delete [] work;
	delete [] rwork;
#endif
	std::cerr << "h_matrix diag: using sunperflib code\n";
    }
#else
    h_matrix H(*this);			     Make a copy of hmx to use
    H.cred(U);				         To Hermitian tridiagonal form
    H.rred(U);				         To real symm. tridiagonal form
    H.tqli(U, D);				     To real diagonal form
    ((n_matrix*)U)->unitary = true;	 This should be unitary now
#endif
    return;				             Now hmx = U*D*adj(U)
}
void h_matrix :: diag ( _matrix * ( & D ) , _matrix * ( & U ) ) { #if defined ( _USING_SUNPERFLIB_ ) || defined ( _USING_LAPACK_ ) int nrows = this -> rows ( ) ; int ncols = this -> cols ( ) ; if ( nrows < 32 ) { h_matrix H ( * this ) ; H . cred ( U ) ; H . rred ( U ) ; H . tqli ( U , D ) ; ( ( n_matrix * ) U ) -> unitary = true ; } else { n_matrix * hmx = new n_matrix ( nrows , ncols ) ; this -> convert ( hmx ) ; n_matrix * hmx1 = ( n_matrix * ) hmx -> transpose ( ) ; #ifdef _USING_SUNPERFLIB_ char jobz = ' ' ; char uplo = ' ' ; int N = nrows ; int lda = nrows ; double * w_eig = new double [ N ] ; int info = -55555 ; zheev ( jobz , uplo , N , ( doublecomplex * ) hmx1 -> data , lda , w_eig , & info ) ; #endif #ifdef _USING_LAPACK_ char jobz = ' ' ; char uplo = ' ' ; #if defined ( __LP64__ ) int N = nrows ; int lda = nrows ; int info = -55555 ; int lwork = 2 * N + 10 ; #else long int N = nrows ; long int lda = nrows ; long int info = -55555 ; long int lwork = 2 * N + 10 ; #endif double * w_eig = new double [ N ] ; complex * work = new complex [ 2 * lwork ] ; double * rwork = new double [ 3 * N - 2 ] ; zheev_ ( & jobz , & uplo , & N , ( __CLPK_doublecomplex * ) hmx1 -> data , & lda , w_eig , ( __CLPK_doublecomplex * ) work , & lwork , rwork , & info ) ; #endif if ( info == 0 ) { } else if ( info > 0 ) { std :: cerr << " \n \n " ; } else if ( info == -55555 ) { std :: cerr << " \n \n " ; } else { std :: cerr << " \n " << info << " \n " ; } D = new d_matrix ( nrows , ncols ) ; U = new n_matrix ( nrows , ncols ) ; complex * hmxp = hmx1 -> data ; for ( int i = 0 ; i < nrows ; i ++ ) { D -> put ( w_eig [ i ] , i , i ) ; for ( int j = 0 ; j < ncols ; j ++ ) { U -> put ( hmxp [ j * nrows + i ] , i , j ) ; } } ( ( n_matrix * ) U ) -> unitary = true ; delete hmx ; delete hmx1 ; delete [ ] w_eig ; #ifdef _USING_LAPACK_ delete [ ] work ; delete [ ] rwork ; #endif } #else h_matrix H ( * this ) ; H . cred ( U ) ; H . rred ( U ) ; H . tqli ( U , D ) ; ( ( n_matrix * ) U ) -> unitary = true ; #endif return ; }
The routine diag combines the routines cred,   rred and tqli (using
the function sign) to diagonalize a complex Hermitian matrix.
The routine diag combines the routines cred rred and tqli ( using the function sign ) to diagonalize a complex Hermitian matrix .
// Input        hmx     : h_matrix (this) //              D	    : Pointer to mx to become diagonal mx //              U    	: Pointer to mx to become eigenvectors mx // // Output       void    : The matrices D and U are set to the diagonal //                          matrix of hmx eigenvalues and the matrix //                          of dmx eigenvectors respectively. // // Note 1               : Both U and D MUST be input pointing to //                          NULL matrix. Their data will be allocated in //                          cred and tqli respectively. Their reference //                          counting MUST be handled externally. // Note 2               : D will be a d_matrix, U n_matrix unitary // then do it the old fashioned gamma way. // Make a copy of hmx to use // To Hermitian tridiagonal form // To real symm. tridiagonal form // To real diagonal form // This should be unitary now //	std::cerr << "h_matrix diag: using gamma code\n"; // Calculate "eigenvectors" AND "eigenvalues". // Upper triangular... // Calculate "eigenvectors" AND "eigenvalues". // Upper triangular... /* this is needed on MacOS CLAPACK to make types match*/ // all is well. // Create the matrixes D and U. // Copy results back into D and U. // This should be unitary now //	std::cerr << "h_matrix diag: using sunperflib code\n"; // Make a copy of hmx to use // To Hermitian tridiagonal form // To real symm. tridiagonal form // To real diagonal form // This should be unitary now // Now hmx = U*D*adj(U)