00001
00002
00003
00004
00005
00006
00007 #ifndef SYNAPS_LINALG_DET_H
00008 #define SYNAPS_LINALG_DET_H
00009
00010 #include <synaps/init.h>
00011 #include <synaps/linalg/MethodName.h>
00012 #include <synaps/linalg/MATRIX.m>
00013 #include <synaps/linalg/Decomp.h>
00014
00015 __BEGIN_NAMESPACE_SYNAPS
00017
00018 template<class X> bool IsNull(const X & x) {return( x==0);}
00019
00020 inline int eps(unsigned int k)
00021 {
00022 if(k%2==0) return 1; else return -1;
00023 }
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 template <class R>
00064 typename R::size_type pseudo_triang( R & M)
00065 {
00066 typedef typename R::size_type size_type;
00067
00068 size_type * colind=new size_type[M.nbcol()], * rowind=new size_type[M.nbrow()];
00069 size_type nbtrans=0;
00070
00071 for (size_type i=0; i<M.nbrow(); i++) {
00072 colind[i]=0;
00073 rowind[i]=0;
00074 }
00075
00076 for (size_type i=0; i<M.nbrow(); i++)
00077 for (size_type j=0; j<M.nbcol(); j++)
00078 if (IsNull(M(i,j))) {
00079 rowind[i]+=1;
00080 colind[j]+=1;
00081 }
00082
00083 bool finish=false;
00084 while (!finish) {
00085 finish=true;
00086 for (size_type i=0; i<M.nbrow()-1; i++)
00087 {
00088 if (rowind[i]>rowind[i+1]) {
00089 M.swaprow(i,i+1);
00090 swap(rowind[i],rowind[i+1]);
00091 nbtrans++;
00092 finish=false;
00093 }
00094 }
00095 }
00096
00097 finish=false;
00098 while (!finish) {
00099 finish=true;
00100 for (size_type i=0; i<M.nbcol()-1; i++)
00101 {
00102 if (colind[i]<colind[i+1]) {
00103 M.swapcol(i,i+1);
00104 swap(colind[i],colind[i+1]);
00105 nbtrans++;
00106 finish=false;
00107 }
00108 }
00109 }
00110 return nbtrans;
00111 }
00112
00124 template <class R>
00125 typename R::value_type Det(const R & M){
00126 using namespace MATRIX;
00127 assert(M.nbcol()==M.nbrow());
00128
00129 if(M.nbcol()==2)
00130 return M(0,0)*M(1,1)-M(1,0)*M(0,1);
00131 else if(M.nbcol()==3) {
00132 typename R::value_type res;
00133 res = M(0,0)*M(1,1)*M(2,2);
00134 res -=M(0,0)*M(1,2)*M(2,1);
00135 res -=M(1,0)*M(0,1)*M(2,2);
00136 res +=M(1,0)*M(0,2)*M(2,1);
00137 res +=M(2,0)*M(0,1)*M(1,2);
00138 res -=M(2,0)*M(0,2)*M(1,1);
00139 return res;
00140 } else
00141 {
00142 R N(M.rep());
00143 return decomp(Bareiss(),N.rep());
00144 }
00145 }
00146
00147
00159 template <class R>
00160 typename R::value_type Det(const Mixed & mth, R & M){
00161 assert(M.nbcol()==M.nbrow());
00162
00163 typedef typename R::value_type value_type;
00164 typedef typename R::size_type size_type;
00165
00166 value_type delta(0);
00167 size_type opp(0);
00168
00169 if (M.nbcol()==1) return M(0,0);
00170 if (M.nbcol()==2)
00171 return (M(0,0)*M(1,1)-M(0,1)*M(1,0));
00172
00173 opp=pseudo_triang(M);
00174
00175 for(size_type k=0; k< M.nbcol();k++) {
00176 if (!IsNull(M(0,k))) {
00177 R N(M.nbrow(),M.nbcol());
00178 for(size_type i=0; i< M.nbrow();i++)
00179 for(size_type j=0; j< M.nbcol();j++)
00180 N(i,j) = M(i,j);
00181
00182 submatrix(N,0,0,k,k);
00183
00184 if ( M.nbcol() >= 4) {
00185 bool prob=true, altern=false;
00186 size_type i=0;
00187 for (size_type j=0;j<N.nbrow();j++) {
00188 while ((prob) && (i<=j)) {
00189 if (IsNull(N(j,i))) i++;
00190 else prob=false;
00191 }
00192 if (prob) {
00193 size_type l=j+1;
00194 while (l<N.nbrow() && IsNull(N(l,j)) ) l++;
00195 if (l<N.nbrow()) N.swaprow(j,l);
00196 else if (((j+1)<N.nbrow()) && (!IsNull(N(j,j+1)))) N.swapcol(j,j+1);
00197 if (!altern) altern=true;
00198 else altern=false;
00199 }
00200 i=0;
00201 prob=true;
00202 }
00203
00204 Triang<Bareiss>(N);
00205 if (altern) delta-= value_type(eps(k))*M(0,k)*N(M.nbrow()-2,M.nbcol()-2);
00206 else delta+= value_type(eps(k))*M(0,k)*N(M.nbrow()-2,M.nbcol()-2);
00207 }
00208 else
00209 if (M.nbcol()==4) {
00210 value_type minor1,minor2,minor3;
00211 minor1 =N(1,1)*N(2,2);
00212 minor1-=N(1,2)*N(2,1);
00213 minor1*=N(0,0);
00214 minor2 =N(1,0)*N(2,2);
00215 minor2-=N(1,2)*N(2,0);
00216 minor2*=N(0,1);
00217 minor2*=(value_type)(-1);
00218 minor3 =N(1,0)*N(1,2);
00219 minor3-=N(1,1)*N(2,0);
00220 minor3*=N(0,2);
00221 delta+=value_type(eps(k))*M(0,k)*(minor1-minor2+minor3);
00222 }
00223 else
00224 if (M.nbcol()==3)
00225 delta+= value_type(eps(k))*M(0,k)*(N(0,0)*N(1,1)-N(0,1)*N(1,0));
00226 }
00227 }
00228
00229 if (opp%2 == 0)
00230 return delta;
00231 else
00232 return -delta;
00233 }
00234
00244 template<class R>
00245 typename R::value_type Det(const Berko & mth, const R & M)
00246 {
00247 using MATRIX::trace;
00248 assert(M.nbcol()==M.nbrow());
00249 typedef typename R::value_type value_type;
00250 typedef typename R::size_type size_type;
00251
00252 R SM(M);
00253 value_type d; int eps=1;
00254 d=trace(SM);
00255 for (size_type i=2; i<=M.nbcol(); i++)
00256 {
00257 SM*=M;
00258 SM-=(M*d);
00259 d=trace(SM);
00260 d/=i; eps*=-1;
00261 }
00262 d*=eps;
00263 return d;
00264 }
00265
00266
00267 __END_NAMESPACE_SYNAPS
00268
00269 #endif // SYNAPS_LINALG_DET_H
00270
00271