5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
12 /*********************************************************************************************************************************/
14 inline double SIGN(const double a, const double b)
16 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
18 /*********************************************************************************************************************************/
20 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
22 vector<vector<double> > product;
24 int first_rows = first.size();
25 int first_cols = first[0].size();
26 int second_cols = second[0].size();
28 product.resize(first_rows);
29 for(int i=0;i<first_rows;i++){
30 product[i].resize(second_cols);
33 for(int i=0;i<first_rows;i++){
34 for(int j=0;j<second_cols;j++){
36 if (m->control_pressed) { return product; }
39 for(int k=0;k<first_cols;k++){
40 product[i][j] += first[i][k] * second[k][j];
48 m->errorOut(e, "LinearAlgebra", "matrix_mult");
54 /*********************************************************************************************************************************/
56 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
58 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
60 double scale, hh, h, g, f;
67 for(int i=n-1;i>0;i--){
71 for(int k=0;k<l+1;k++){
72 scale += fabs(a[i][k]);
78 for(int k=0;k<l+1;k++){
80 h += a[i][k] * a[i][k];
83 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
88 for(int j=0;j<l+1;j++){
89 a[j][i] = a[i][j] / h;
91 for(int k=0;k<j+1;k++){
92 g += a[j][k] * a[i][k];
94 for(int k=j+1;k<l+1;k++){
95 g += a[k][j] * a[i][k];
101 for(int j=0;j<l+1;j++){
103 e[j] = g = e[j] - hh * f;
104 for(int k=0;k<j+1;k++){
105 a[j][k] -= (f * e[k] + g * a[i][k]);
120 for(int i=0;i<n;i++){
123 for(int j=0;j<l;j++){
125 for(int k=0;k<l;k++){
126 g += a[i][k] * a[k][j];
128 for(int k=0;k<l;k++){
129 a[k][j] -= g * a[k][i];
135 for(int j=0;j<l;j++){
136 a[j][i] = a[i][j] = 0.0;
142 catch(exception& e) {
143 m->errorOut(e, "LinearAlgebra", "tred2");
148 /*********************************************************************************************************************************/
150 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
152 /*********************************************************************************************************************************/
154 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
156 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
159 double s, r, p, g, f, dd, c, b;
162 for(int i=1;i<=n;i++){
167 for(int l=0;l<n;l++){
170 for(myM=l;myM<n-1;myM++){
171 dd = fabs(d[myM]) + fabs(d[myM+1]);
172 if(fabs(e[myM])+dd == dd) break;
175 if(iter++ == 30) cerr << "Too many iterations in tqli\n";
176 g = (d[l+1]-d[l]) / (2.0 * e[l]);
178 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
181 for(i=myM-1;i>=l;i--){
184 e[i+1] = (r=pythag(f,g));
193 r = (d[i] - g) * s + 2.0 * c * b;
194 d[i+1] = g + ( p = s * r);
196 for(int k=0;k<n;k++){
198 z[k][i+1] = s * z[k][i] + c * f;
199 z[k][i] = c * z[k][i] - s * f;
202 if(r == 0.00 && i >= l) continue;
211 for(int i=0;i<n;i++){
213 for(int j=i;j<n;j++){
221 for(int j=0;j<n;j++){
231 catch(exception& e) {
232 m->errorOut(e, "LinearAlgebra", "qtli");
236 /*********************************************************************************************************************************/