src/mmdb_math_linalg.h
Go to the documentation of this file.
00001 // $Id: mmdb_math_linalg.h $ 00002 // ================================================================= 00003 // 00004 // CCP4 Coordinate Library: support of coordinate-related 00005 // functionality in protein crystallography applications. 00006 // 00007 // Copyright (C) Eugene Krissinel 2000-2013. 00008 // 00009 // This library is free software: you can redistribute it and/or 00010 // modify it under the terms of the GNU Lesser General Public 00011 // License version 3, modified in accordance with the provisions 00012 // of the license to address the requirements of UK law. 00013 // 00014 // You should have received a copy of the modified GNU Lesser 00015 // General Public License along with this library. If not, copies 00016 // may be downloaded from http://www.ccp4.ac.uk/ccp4license.php 00017 // 00018 // This program is distributed in the hope that it will be useful, 00019 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00021 // GNU Lesser General Public License for more details. 00022 // 00023 // ================================================================= 00024 // 00025 // 12.09.13 <-- Date of Last Modification. 00026 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 00027 // ----------------------------------------------------------------- 00028 // 00029 // **** Module : LinAlg <interface> 00030 // ~~~~~~~~~ 00031 // **** Project : MMDB ( MacroMolecular Data Base ) 00032 // ~~~~~~~~~ 00033 // 00034 // (C) E.Krissinel 2000-2013 00035 // 00036 // ================================================================= 00037 // 00038 // 00039 00040 #ifndef __MMDB_MATH_LinAlg__ 00041 #define __MMDB_MATH_LinAlg__ 00042 00043 #include "mmdb_mattype.h" 00044 00045 namespace mmdb { 00046 00047 namespace math { 00048 00049 // ========================== Jacobi ============================= 00050 00051 00054 extern void Jacobi ( int N, 00055 rmatrix A, 00056 00057 00058 rmatrix T, 00059 rvector Eigen, 00060 00061 rvector Aik, 00062 int & Signal 00063 00064 ); 00065 00066 00067 // A5.5.2 : Perturbed Cholesky Decomposition 00068 extern void PbCholDecomp ( int N, 00069 rvector HDiag, 00070 realtype MaxOff, 00071 realtype MachEps, 00072 rmatrix L, 00073 realtype & MaxAdd ); 00074 00075 // A3.2.3a : Cholesky's L - Solution of 00076 // L*Y = B ( given B ) 00077 extern void LSolve ( int N, rmatrix L, rvector B, rvector Y ); 00078 00079 // A3.2.3b : Cholesky's LT - Solution of 00080 // LT*X = Y ( given Y ) 00081 extern void LTSolve ( int N, rmatrix L, rvector Y, rvector X ); 00082 00083 // A3.2.3 : Solution of the equation L*LT*S = G 00084 // by the Cholesky's method 00085 extern void ChSolve ( int N, rmatrix L, rvector G, rvector S ); 00086 00087 00088 // ---------------------------------------------------- 00089 00090 extern void FastInverse ( int N, rmatrix A, ivector J0, 00091 //#D realtype & Det, 00092 int & Signal ); 00093 // 00094 // 13.09.90 <-- Last Modification Date 00095 // ------------------------ 00096 // 00097 // ================================================ 00098 // 00099 // Fast Inversion of the matrix A 00100 // by the method of GAUSS - JORDAN . 00101 // 00102 // ------------------------------------------------ 00103 // 00104 // Input parameters are : 00105 // 00106 // N - dimension of the matrix 00107 // A - the matrix [1..N][1..N] to be inverted. 00108 // ------------------------------------------------ 00109 // 00110 // J0 - integer vector [1..N] for temporal storage 00111 // 00112 // ------------------------------------------------ 00113 // 00114 // Output parameters are : 00115 // 00116 // A - the inverted matrix 00117 // Signal - the error key : 00118 // = 0 <=> O'K 00119 // else 00120 // degeneration was found, and 00121 // the rang of matrix is Signal-1. 00122 // 00123 // Variable Det may return the determinant 00124 // of matrix A. To obtain it, remove all comments 00125 // of form //#D. 00126 // 00127 // ================================================ 00128 00129 00130 // ---------------------------------------------------- 00131 00132 extern void SVD ( int NA, int M, int N, 00133 rmatrix A, rmatrix U, rmatrix V, 00134 rvector W, rvector RV1, 00135 bool MatU, bool MatV, 00136 int & RetCode ); 00137 // 00138 // 13.12.01 <-- Last Modification Date 00139 // ------------------------ 00140 // 00141 // ================================================ 00142 // 00143 // The Singular Value Decomposition 00144 // of the matrix A by the algorithm from 00145 // G.Forsait, M.Malkolm, K.Mouler. Numerical 00146 // methods of mathematical calculations // 00147 // M., Mir, 1980. 00148 // 00149 // Matrix A is represented as 00150 // 00151 // A = U * W * VT 00152 // 00153 // ------------------------------------------------ 00154 // 00155 // All dimensions are indexed from 1 on. 00156 // 00157 // ------------------------------------------------ 00158 // 00159 // Input parameters: 00160 // 00161 // NA - number of lines in A. NA may be 00162 // equal to M or N only. If NA=M 00163 // then usual SVD will be made. If MA=N 00164 // then matrix A is transposed before 00165 // the decomposition, and the meaning of 00166 // output parameters U and V is 00167 // swapped (U accepts VT and VT accepts U). 00168 // In other words, matrix A has physical 00169 // dimension of M x N , same as U and V; 00170 // however the logical dimension of it 00171 // remains that of N x M . 00172 // M - number of lines in U 00173 // N - number of columns in U,V and length 00174 // of W,RV1 . Always provide M >= N ! 00175 // A - matrix [1..M][1..N] or [1..N][1..M] 00176 // to be decomposed. The matrix does not 00177 // change, and it may coincide with U or 00178 // V, if NA=M (in which case A does change) 00179 // MatU - compute U , if set True 00180 // MatV - compute V , if set True 00181 // RV1 - temporary array [1..N]. 00182 // U - should be always supplied as an array of 00183 // [1..M][1..N], M>=N . 00184 // V - should be suuplied as an array of 00185 // [1..N][1..N] if MatV is True . 00186 // 00187 // ------------------------------------------------ 00188 // 00189 // Output parameters are : 00190 // 00191 // W - N non-ordered singular values, 00192 // if RetCode=0. If RetCode<>0, the 00193 // RetCode+1 ... N -th values are still 00194 // valid 00195 // U - matrix of right singular vectors 00196 // (arranged in columns), corresponding 00197 // to the singular values in W, if 00198 // RetCode=0 and MatU is True. If MatU 00199 // is False, U is still used as a 00200 // temporary array. If RetCode<>0 then 00201 // the RetCode+1 ... N -th vectors 00202 // are valid 00203 // V - matrix of left singular vectors 00204 // (arranged in columns), corresponding 00205 // to the singular values in W, if 00206 // RetCode=0 and MatV is True. If MatV 00207 // is False, V is not used and may be set 00208 // to NULL. If RetCode<>0 then the 00209 // RetCode+1 ... N -th vectors are valid 00210 // RetCode - the error key : 00211 // = 0 <=> O'K 00212 // else 00213 // = k, if the k-th singular value 00214 // was not computed after 30 iterations. 00215 // 00216 // ------------------------------------------------ 00217 // 00218 // Key Variables are : 00219 // 00220 // ItnLimit - the limit for iterations 00221 // 00222 // This routine does not use any machine-dependent 00223 // constants. 00224 // 00225 // ================================================ 00226 // 00227 // 00228 00229 extern void OrderSVD ( int M, int N, rmatrix U, rmatrix V, 00230 rvector W, bool MatU, bool MatV ); 00231 00232 00233 } 00234 } 00235 00236 #endif
