00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef TRISLV_H
00029 #define TRISLV_H
00030
00031
00032 #include "tnt_triang.h"
00033
00034 namespace SPUC
00035 {
00036
00037 template <class MaTriX, class VecToR>
00038 VecToR lower_triangular_solve( MaTriX &A, VecToR &b)
00039 {
00040 int N = A.dim1();
00041
00042
00043
00044 assert(A.dim2() == N);
00045 assert(b.dim() == N);
00046
00047 VecToR x(N);
00048
00049 int i;
00050 for (i=1; i<=N; i++)
00051 {
00052 typename MaTriX::element_type tmp=0;
00053
00054 for (int j=1; j<i; j++)
00055 tmp = tmp + A(i,j)*x(j);
00056
00057 x(i) = (b(i) - tmp)/ A(i,i);
00058 }
00059
00060 return x;
00061 }
00062
00063
00064 template <class MaTriX, class VecToR>
00065 VecToR unit_lower_triangular_solve( MaTriX &A, VecToR &b)
00066 {
00067 int N = A.dim1();
00068
00069
00070
00071 assert(A.dim2() == N);
00072 assert(b.dim() == N);
00073
00074 VecToR x(N);
00075
00076 int i;
00077 for (i=1; i<=N; i++)
00078 {
00079
00080 typename MaTriX::element_type tmp=0;
00081
00082 for (int j=1; j<i; j++)
00083 tmp = tmp + A(i,j)*x(j);
00084
00085 x(i) = b(i) - tmp;
00086 }
00087
00088 return x;
00089 }
00090
00091
00092 template <class MaTriX, class VecToR>
00093 VecToR linear_solve( lowerTriangularView<MaTriX> &A,
00094 VecToR &b)
00095 {
00096 return lower_triangular_solve(A, b);
00097 }
00098
00099 template <class MaTriX, class VecToR>
00100 VecToR linear_solve( unitlowerTriangularView<MaTriX> &A,
00101 VecToR &b)
00102 {
00103 return unit_lower_triangular_solve(A, b);
00104 }
00105
00106
00107
00108
00109
00110 template <class MaTriX, class VecToR>
00111 VecToR upper_triangular_solve( MaTriX &A, VecToR &b)
00112 {
00113 int N = A.dim1();
00114
00115
00116
00117 assert(A.dim2() == N);
00118 assert(b.dim() == N);
00119
00120 VecToR x(N);
00121
00122 int i;
00123 for (i=N; i>=1; i--)
00124 {
00125
00126 typename MaTriX::element_type tmp=0;
00127
00128 for (int j=i+1; j<=N; j++)
00129 tmp = tmp + A(i,j)*x(j);
00130
00131 x(i) = (b(i) - tmp)/ A(i,i);
00132 }
00133
00134 return x;
00135 }
00136
00137
00138 template <class MaTriX, class VecToR>
00139 VecToR unit_upper_triangular_solve( MaTriX &A, VecToR &b)
00140 {
00141 int N = A.dim1();
00142
00143
00144
00145 assert(A.dim2() == N);
00146 assert(b.dim() == N);
00147
00148 VecToR x(N);
00149
00150 int i;
00151 for (i=N; i>=1; i--)
00152 {
00153
00154 typename MaTriX::element_type tmp=0;
00155
00156 for (int j=i+1; j<i; j++)
00157 tmp = tmp + A(i,j)*x(j);
00158
00159 x(i) = b(i) - tmp;
00160 }
00161
00162 return x;
00163 }
00164
00165
00166 template <class MaTriX, class VecToR>
00167 VecToR linear_solve( upperTriangularView<MaTriX> &A,
00168 VecToR &b)
00169 {
00170 return upper_triangular_solve(A, b);
00171 }
00172
00173 template <class MaTriX, class VecToR>
00174 VecToR linear_solve( unitupperTriangularView<MaTriX> &A,
00175 VecToR &b)
00176 {
00177 return unit_upper_triangular_solve(A, b);
00178 }
00179
00180
00181 }
00182
00183 #endif
00184