$search
00001 #ifndef MATRICES_H 00002 #define MATRICES_H 00003 #include <math.h> 00004 00005 typedef float** matrix; // designate a matrix as an array (2D) of floats 00006 00007 void addMatrix (matrix A, matrix B, matrix C, int m, int n); 00008 void addVector (float* A, float* B, float* C, int m); 00009 void subMatrix (matrix A, matrix B, matrix C, int m, int n); 00010 void subVector (float* A, float* B, float* C, int m); 00011 void negMatrix(matrix A, int m, int n); 00012 void multMatrix(matrix A, matrix B, matrix C, int Am, int An, int Bn); 00013 void multMatVec(matrix A, float* B, float* C, int Am, int Bn); 00014 void transpose(matrix A, matrix B, int Am, int An); 00015 void invert3(matrix A, matrix B); 00016 void invert2(matrix A, matrix B); 00017 float determinant(matrix A,matrix tmp, int n); 00018 void getMinor(matrix A, matrix B,int An, int m, int n); 00019 void invt(matrix A, matrix B, matrix tmp1, matrix tmp2, int An); 00020 float BhQhBhT(matrix Qh, float* Bh, float* tmp); 00021 void VectdivSc(float* V, float scl); 00022 00023 00024 void addMatrix (matrix A, matrix B, matrix C, int m, int n){ 00025 /*C= A+B; 00026 * A,B and C need to be the same size (m by n); function does not check*/ 00027 int i, j; 00028 for(i=0; i<m; ++i){ 00029 for(j=0; j<n; ++j){ 00030 C[i][j] = A[i][j] + B[i][j]; 00031 }//inner for 00032 }//outer for 00033 }//add matrices 00034 00035 void addVector (float* A, float* B, float* C, int m){ 00036 /*C= A+B; 00037 * A,B and C need to be the same size (m); function does not check*/ 00038 int i; 00039 for(i=0; i<m; ++i){ 00040 C[i] = A[i] + B[i]; 00041 }//for 00042 }//add Vectors 00043 00044 void subMatrix (matrix A, matrix B, matrix C, int m, int n){ 00045 /*C= A-B; 00046 * A,B and C need to be the same size (m by n); function does not check*/ 00047 int i, j; 00048 for(i=0; i<m; ++i){ 00049 for(j=0; j<n; ++j){ 00050 C[i][j] = A[i][j] - B[i][j]; 00051 }//inner for 00052 }//outer for 00053 }//sub matrices 00054 00055 void subVector (float* A, float* B, float* C, int m){ 00056 /*C= A-B; 00057 * A,B and C need to be the same size (m); function does not check*/ 00058 int i; 00059 for(i=0; i<m; ++i){ 00060 C[i] = A[i] - B[i]; 00061 }//for 00062 }//add Vectors 00063 00064 00065 00066 void negMatrix(matrix A, int m, int n){ 00067 /*negate matrix*/ 00068 int i,j; 00069 for(i=0; i<m; ++i){ 00070 for(j=0; j<n; ++j){ 00071 A[i][j] = -A[i][j]; 00072 }//inner for 00073 }//outer for 00074 }//negMAtrix 00075 00076 00077 void multMatrix(matrix A, matrix B, matrix C, int Am, int An, int Bn){ 00078 //C = A*B 00079 //A,B and C need to be of correct dimensions, function does not check 00080 int i,j,k; 00081 float sum; 00082 00083 for(i=0; i < Bn; i++){ 00084 for(j=0; j < Am; j++){ 00085 sum=0; 00086 for(k=0; k < An; k++) 00087 sum += A[j][k] * B[k][i]; 00088 C[j][i] = sum; 00089 }//inner for 00090 }//outer for 00091 } 00092 00093 void multMatVec(matrix A, float* B, float* C, int Am, int Bn){ 00094 /*C = A*B 00095 *A,B and C need to be of correct dimensions, function does not check */ 00096 00097 00098 int j,k; 00099 float sum; 00100 for(j=0; j < Am; j++){ 00101 sum=0; 00102 for(k=0; k < Bn; k++) 00103 sum += A[j][k] * B[k]; 00104 C[j] = sum; 00105 }//for 00106 }//multMatVec 00107 00108 00109 void transpose(matrix A, matrix B, int Am, int An){ 00110 /*B = A'*/ 00111 int i,j; 00112 for(i=0; i < Am; i++) 00113 for(j=0; j < An; j++) 00114 B[j][i] = A[i][j]; 00115 00116 } 00117 00118 void invert3(matrix A, matrix B){ 00119 /* B = A^-1 for 3x3 matrices*/ 00120 float det = (A[0][0]*A[1][1]*A[2][2]-A[0][0]*A[1][2]*A[2][1]-A[1][0]*A[0][1]*A[2][2]+A[1][0]*A[0][2]*A[2][1]+A[2][0]*A[0][1]*A[1][2]-A[2][0]*A[0][2]*A[1][1]); 00121 00122 B[0][0]=(A[1][1]*A[2][2]-A[1][2]*A[2][1])/det; 00123 B[0][1]=(-A[0][1]*A[2][2]+A[0][2]*A[2][1])/det; 00124 B[0][2]=(A[0][1]*A[1][2]-A[0][2]*A[1][1])/det; 00125 B[1][0]=(-A[1][0]*A[2][2]+A[1][2]*A[2][0])/det; 00126 B[1][1]=(A[0][0]*A[2][2]-A[0][2]*A[2][0])/det; 00127 B[1][2]=(-A[0][0]*A[1][2]+A[0][2]*A[1][0])/det; 00128 B[2][0]=(A[1][0]*A[2][1]-A[1][1]*A[2][0])/det; 00129 B[2][1]=(-A[0][0]*A[2][1]+A[0][1]*A[2][0])/det; 00130 B[2][2]=(A[0][0]*A[1][1]-A[0][1]*A[1][0])/det; 00131 }//invert3 00132 00133 void invert2(matrix A, matrix B){ 00134 /* B = A^-1 for 2x2 matrices*/ 00135 float det = (A[0][0]*A[1][1]-A[0][1]*A[1][0]); 00136 00137 B[0][0]=(A[1][1])/det; 00138 B[0][1]=(-A[0][1])/det; 00139 B[1][0]=(-A[1][0])/det; 00140 B[1][1]=(A[0][0])/det; 00141 }//invert2 00142 00143 float determinant(matrix A,matrix tmp, int n){ 00144 /*returns determinant of n x n matrix*/ 00145 int i,j,k; 00146 float q; 00147 float dtr=1; 00148 /* copy A to tmp*/ 00149 for(i=0; i < n; i++) 00150 for(j=0; j < n; j++) 00151 tmp[i][j] = A[j][i]; 00152 00153 /*compute determinant*/ 00154 for(i=0;i < n;i++){ 00155 for(j=0;j < n;j++){ 00156 q=tmp[j][i]/tmp[i][i]; 00157 for(k=0;k < n;k++){ 00158 if(i==j) break; 00159 tmp[j][k]=(tmp[j][k]-tmp[i][k])*q; 00160 }//innermost for 00161 }//middle for 00162 }//outer for 00163 for(i=0;i < n;i++){ 00164 dtr=dtr*tmp[i][i]; 00165 }//for 00166 return(dtr); 00167 }//det 00168 00169 00170 void getMinor(matrix A, matrix B,int An, int m, int n){ 00171 /* use this function to generate a matrix B that is the mn minor of A; A must be square 00172 * m>=1, n>=1, An is the size of A*/ 00173 int i,j; 00174 int k=0; 00175 for(i=0; i < An; i++){ 00176 if(i!=(m-1)){ 00177 for(j=0; j < (An-1); j++){ 00178 if(j<(n-1)){ 00179 B[k][j] = A[i][j]; 00180 } 00181 else{ 00182 B[k][j] = A[i][j+1]; 00183 } 00184 }//inner for (column) 00185 k++; 00186 }//if 00187 }//outer for 00188 }//get minor 00189 00190 00191 void invt(matrix A, matrix B, matrix tmp1, matrix tmp2, int An){ 00192 /*B = A^-1 ; tmp1 & tmp2 must be at least of size (An-1)x(An-1) 00193 * function does not check if matrix is singular, tried, but even if det(A)==0, it did not escape, 00194 * just reutrned "nan" for some elements, and garbage for others 00195 * */ 00196 int i,j; 00197 float detA = determinant(A,tmp2,An); 00198 for(i=0; i < An; i++){ 00199 for(j=0; j < An; j++){ 00200 getMinor(A,tmp1,An,i+1,j+1); 00201 if(((i+j)%2)==0){ 00202 B[j][i]= determinant(tmp1,tmp2,(An-1))/detA; 00203 } 00204 else{ 00205 B[j][i]= -determinant(tmp1,tmp2,(An-1))/detA; 00206 } 00207 //problem if an element is zero (returns "nan" for zero element) 00208 }//inner for 00209 }//outer for 00210 }//invert 00211 00212 float BhQhBhT(matrix Qh, float* Bh, float* tmp){ 00213 //function returns the value for Bh*Qh*Bh', where Qh is an 3x3 matrix, and Bh is a vector of length 3 00214 00215 00216 00217 //Bh*Qh 00218 int i,j; 00219 float sum; 00220 for(i=0; i < 3; i++){ 00221 sum=0; 00222 for(j=0; j < 3; j++) 00223 sum += Bh[j]*Qh[j][i]; 00224 tmp[i] = sum; 00225 }//outer for 00226 00227 //tmp*BhT 00228 sum = 0; 00229 for(i=0; i < 3; i++) 00230 sum += tmp[i]*Bh[i]; 00231 00232 00233 00234 return sum; 00235 } 00236 00237 void VectdivSc(float* V, float scl){ 00238 //function returns the value of Vector/scalar where vector is 3x1 00239 00240 int i; 00241 for(i=0; i < 3; i++){ 00242 V[i] = V[i]/scl; 00243 }// for 00244 } 00245 00246 #endif /*MATRICES_H_*/ 00247 00248