00001 /******************************************************************************
00002
00003 Copyright 2008 Departamento de Realidad Virtual
00004 y Unidad de Cómputo Aplicado DGSGA, UNAM.
00005
00006
00007 This file is part of RBF++.
00008
00009 RBF++ is free software: you can redistribute it and/or modify
00010 it under the terms of the GNU General Public License as published by
00011 the Free Software Foundation, either version 3 of the License, or
00012 (at your option) any later version.
00013
00014 RBF++ is distributed in the hope that it will be useful,
00015 but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00017 GNU General Public License for more details.
00018
00019 You should have received a copy of the GNU General Public License
00020 along with RBF++. If not, see <http://www.gnu.org/licenses/>.
00021
00022
00023 *******************************************************************************/
00024
00025 /*
00026 Author: Daniel Cervantes Cabrera.
00027 Project: RBF++
00028 Institution: DGSCA, UNAM.
00029 Date: 3/05/08
00030 Description: Template matrix class.
00031 */
00032
00033
00034 #include <cmath>
00035 #include <cstdlib>
00036 #include <LA/Matrix.h>
00037
00038
00039
00040
00041 namespace LA{
00042
00043 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
00044
00045 template<typename T>
00046 Matrix<T>::Matrix() : data(0), row(0), col(0) {
00047 }
00048
00049 template<typename T>
00050 Matrix<T>::Matrix(int r, int c) :row(r), col(c) {
00051
00052 data = new T* [row];
00053 for(int i = 0; i < row; ++i)
00054 data[i] = new T [col];
00055 }
00056
00057 template<typename T>
00058 Matrix<T>::Matrix(const Matrix<T>& otra) {
00059
00060 if(this != &otra) {
00061 row = otra.row;
00062 col = otra.col;
00063 data = new T* [row];
00064 for(int i = 0; i < row; ++i)
00065 data[i] = new T [col];
00066
00067 for(int i = 0; i < row; ++i)
00068 for(int j = 0; j < col; ++j)
00069 data[i][j] = otra.data[i][j];
00070 }
00071 }
00072
00073 template<typename T>
00074 Matrix<T>::~Matrix() {
00075 for(int i = 0; i < row; ++i)
00076 delete [] data[i];
00077 delete [] data;
00078 }
00079
00080 template<typename T>
00081 void Matrix<T>::resize(int r, int c) {
00082
00083 if( row != 0 && col != 0 ) {
00084 for(int i = 0; i < row; ++i)
00085 delete [] data[i];
00086 delete [] data;
00087 }
00088
00089 row = r;
00090 col = c;
00091
00092 data = new T* [row];
00093 for(int i = 0; i < row; ++i)
00094 data[i] = new T [col];
00095 }
00096
00097 template<typename T>
00098 T& Matrix<T>::operator()(int i, int j) const {
00099 return data[i][j];
00100 }
00101
00102 template<typename T>
00103 Matrix<T> Matrix<T>::transpose() {
00104
00105
00106 Matrix<T> TM(col, row);
00107 for(int i = 0; i < col; ++i)
00108 for(int j = 0; j < row; ++j)
00109 TM(i,j) = data[j][i];
00110
00111 return TM;
00112 }
00113
00114 //Gauss con pivoteo.
00115 template<typename T>
00116 vector<T> Matrix<T>::operator/(vector<T> &b) {
00117
00118 Matrix<T> a;
00119 vector<T> x;
00120 int n = row;
00121 a.resize(n,n+1);
00122 x.resize(n);
00123
00124 T pivot, temp,q;
00125 int ipivot,ip1;
00126
00127 for(int i=0; i<n; i++)
00128 for(int j=0; j<n; j++)
00129 a(i,j) = data[i][j];
00130
00131 for(int i=0; i<n; i++)
00132 a(i,n) = b[i];
00133
00134
00135 for(int i=0; i<n-1; i++){
00136 pivot = 0.0;
00137 for(int j=i; j<n; j++){
00138 temp = fabs(a(j,i));
00139 if(pivot < temp){
00140 pivot = temp;
00141 ipivot = j;
00142 }
00143 }
00144
00145
00146 if(fabs(pivot) < 1e-10)
00147 {
00148 std::cout << "Error matriz singular" << std::endl;
00149 exit(1);
00150 }
00151
00152 if(ipivot != i)
00153 for(int k=i; k <= n; k++)
00154 SWAP(a(i,k),a(ipivot,k));
00155
00156 ip1 = i+1;
00157 for(int k=ip1; k<n; k++){
00158 q = -a(k,i)/a(i,i);
00159 a(k,i) = 0.0;
00160 for(int j = ip1; j<=n; j++)
00161 a(k,j) = q*a(i,j) + a(k,j);
00162 }
00163 }
00164
00165 x[n-1] = a(n-1,n)/a(n-1,n-1);
00166
00167 for(int k=1; k <= n-1; k++)
00168 {
00169 q = 0.0;
00170 for(int j=0; j < k; j++)
00171 q = q + a(n-1-k,n-1-j)*x[n-1-j];
00172
00173 x[n-1-k] = (a(n-1-k,n)-q)/(a(n-1-k,n-1-k));
00174 }
00175
00176 return x;
00177
00178
00179 }
00180
00181 template<typename U>
00182 ostream& operator<<(ostream& s, const Matrix<U>& m) {
00183 for(int i = 0; i <m.row; ++i) {
00184 for(int j = 0; j < m.col; ++j)
00185 s << m(i,j) << "\t";
00186 cout << endl;
00187 }
00188 return s;
00189 }
00190 }
Superior

Tenemos 9 visitantes y ningun miembro en Línea