1 /*************************************************************************** 2 * Copyright (C) 2009 by BUI Quang Minh * 3 * minh.bui@univie.ac.at * 4 * * 5 * This program is free software; you can redistribute it and/or modify * 6 * it under the terms of the GNU General Public License as published by * 7 * the Free Software Foundation; either version 2 of the License, or * 8 * (at your option) any later version. * 9 * * 10 * This program is distributed in the hope that it will be useful, * 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 13 * GNU General Public License for more details. * 14 * * 15 * You should have received a copy of the GNU General Public License * 16 * along with this program; if not, write to the * 17 * Free Software Foundation, Inc., * 18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * 19 ***************************************************************************/ 20 #ifndef EIGENDECOMPOSITION_H 21 #define EIGENDECOMPOSITION_H 22 23 //const double ZERO_FREQ = 0.000001; 24 const double ZERO_FREQ = 1e-10; 25 26 27 /** 28 Eigenvalues, eigenvectors decomposition 29 30 @author BUI Quang Minh <minh.bui@univie.ac.at> 31 */ 32 class EigenDecomposition{ 33 public: 34 EigenDecomposition(); 35 36 ~EigenDecomposition(); 37 38 /** 39 EigenSystem for symmetric matrix 40 @param rate_params rate parameters (not the rate matrix) 41 @param state_freq state frequencies 42 @param eval (OUT) eigenvalues 43 @param evec (OUT) eigenvectors 44 @param inv_evec (OUT) inverse matrix of eigenvectors 45 @param num_state (IN) number of states 46 */ 47 void eigensystem_sym(double **rate_params, double *state_freq, 48 double *eval, double *evec, double *inv_evec, int num_state); 49 50 /** 51 EigenSystem for general non-symmetric matrix 52 @param rate_params rate parameters (not the rate matrix) 53 @param state_freq state frequencies 54 @param eval (OUT) eigenvalues 55 @param evec (OUT) eigenvectors 56 @param inv_evec (OUT) inverse matrix of eigenvectors 57 @param num_state (IN) number of states 58 */ 59 void eigensystem(double **rate_params, double *state_freq, 60 double *eval, double **evec, double **inv_evec, int num_state); 61 62 /** 63 EigenSystem for general non-symmetric matrix without state frequencies 64 @param rate_matrix rate matrix 65 @param eval (OUT) real part of eigenvalues 66 @param eval_imag (OUT) imaginary part of eigenvalues 67 @param evec (OUT) eigenvectors 68 @param inv_evec (OUT) inverse matrix of eigenvectors 69 @param num_state (IN) number of states 70 */ 71 void eigensystem_nonrev(double *rate_matrix, double *state_freq, double *eval, double *eval_imag, 72 double *evec, double *inv_evec, int num_state); 73 74 75 /** TRUE to normalize rate matrix to 1.0 subst per unit time */ 76 bool normalize_matrix; 77 78 /** 79 the total number of substitutions per unit time 80 */ 81 double total_num_subst; 82 83 /** TRUE to ignore state_freq in computation, default: FALSE */ 84 bool ignore_state_freq; 85 86 87 protected: 88 89 /** 90 compute the rate matrix and then normalize it such that the total number of substitutions is 1. 91 @param rate_matrix (IN/OUT) As input, it contains rate parameters. On output it is filled with rate matrix entries 92 @param state_freq state frequencies 93 @param num_state number of states 94 */ 95 virtual void computeRateMatrix(double **rate_matrix, double *state_freq, int num_state); 96 97 /** 98 Eliminate zero entries in the rate matrix. 99 Return the new non-zero matrix with possibly reduced dimension. 100 @param mat input rate matrix 101 @param forg state frequencies 102 @param num number of states 103 @param new_mat (OUT) the new rate matrix 104 @param new_forg (OUT) new state frequencies 105 @param new_num (OUT) new number of states 106 */ 107 void eliminateZero(double **mat, double *forg, int num, 108 double **new_mat, double *new_forg, int &new_num); 109 110 /********************************************************* 111 * aided function for symmetric matrix 112 *********************************************************/ 113 114 /** 115 transform the rate matrix into symmetric form, used for subsequent eigen decomposition 116 @param a (IN/OUT) rate matrix 117 @param stateFrq state frequencies 118 @param stateFrq_sqrt square root of state frequencies 119 @param num_state number of states 120 */ 121 void symmetrizeRateMatrix(double **a, double *stateFrq, double *stateFrq_sqrt, int num_state); 122 123 124 /** 125 Householder transformation of symmetric matrix A 126 to tridiagonal form 127 @param a the input matrix, must be symmetric. On output, 128 a is replaced by the orthogonal matrix effecting the transformation 129 @param n the size of matrix a 130 @param d [0..n-1] returned the diagonal elements of the tridiagonal matrix 131 @param e [0..n-1] returned the off-diagonal elements with e[0]=0 132 */ 133 void tred2(double **a, int n, double *d, double *e); 134 135 /** 136 QL algorithm with implicit shifts to determine eigenvalues and 137 eigenvectors of a real tridiagonal symmetric matrix. 138 @param d [0..n-1] diagonal elements of the tridiagonal matrix. 139 On output d return the eigenvalues. 140 @param e [0..n-1] off-diagonal elements of the tridiagonal matrix, e[0] arbitrary. 141 On output e is destroyed. 142 @param n matrix size 143 @param z must be input as the matrix returned by tred2 144 z[k] return the normalized eigenvector corresponding to d[k] 145 */ 146 void tqli(double *d, double *e, int n, double **z); 147 148 /********************************************************* 149 * aided function for non-symmetric matrix 150 *********************************************************/ 151 152 /** 153 convert a non-symmetric matrix into Hessenberg form with zeros everywhere 154 below the diagonal except for the first sub-diagonal row 155 @param a (IN-OUT) the matrix 156 @param ordr (OUT) the order of columns 157 @param n (IN) size of matrix 158 */ 159 void elmhes(double **a, int *ordr, int n); 160 161 /* 162 something here 163 */ 164 void eltran(double **a, double **zz, int *ordr, int n); 165 166 /* 167 something here 168 */ 169 void mcdiv(double ar, double ai, double br, double bi, 170 double *cr, double *ci); 171 172 /** 173 QR algorithm for non-symmetric matrix to calculate eigenvectors and eigenvalues 174 of a Hessenberg matrix (should be preceded by elmhes function) 175 @param n (IN) size of matrix 176 */ 177 void hqr2(int n, int low, int hgh, double **h, double **zz, double *wr, double *wi); 178 179 /** 180 compute the inverse of a square matrix 181 @param inmat (IN) the matrix 182 @param imtrx (OUT) the inverse of the input matrix 183 @param size the size of matrix 184 */ 185 void luinverse(double **inmat, double **imtrx, int size); 186 187 void checkevector(double *evec, double *ivec, int nn); 188 189 }; 190 191 #endif 192