matrix2.txt

几个算法程序

int p, q, count; 
double          a, b, c, d; 
double sintheta, costheta, sin2theta, cos2theta; 
double *bb=NULL; // temporary buffer 
double *tip, *tiq, *tpj, *tqj; 
double tpp, tpq, tqp, tqq; 
 
bb=new double[nn*4]; 
if ( bb==NULL ) 
return; 
tip=&bb[0]; 
tiq=&bb[nn]; 
tpj=&bb[2*nn]; 
tqj=&bb[3*nn]; 
 
// S=I 
n=nn*nn; 
for ( i=0;  i<n;  i++ ) 
ss[i]=0.0; 
for ( i=0;  i<nn;  i++ ) 
ss[i*nn+i]=1.0; 
 
count=n; 
while ( --count>0 ) 
{ 
// find max{ |aa[p][q]| } , p!=q 
a=0.0; 
p=0; 
q=0; 
for ( i=0;  i<nn-1;  i++ ) 
{ 
k = i*nn; 
for ( j=i+1;  j<nn;  j++ ) 
{ 
l = k + j; 
if ( fabs(aa[l])>fabs(a) ) 
{ 
a=aa[l]; 
p=i; 
q=j; 
} 
} 
} 
if ( p==q ) 
break; 
// |aa[p][q]|<epsilon: finish 
if ( fabs(a)<MATRIXEPSILON ) 
break; 
 
// calculate angle rotated 
b=-aa[p*nn+q]; 
c=0.5*(aa[q*nn+q] - aa[p*nn+p]); 
d=sqrt(b*b+c*c); 
sin2theta=b/d; 
cos2theta=c/d; 
sintheta=sqrt( (1.0-cos2theta)*0.5 ); 
costheta=sqrt( (cos2theta+1.0)*0.5 ); 
// 0<=2theta<2pi, 0<=theta<pi 
// if ( b<0.0 )// sin2theta<0.0: pi<=2theta<2pi, pi/2<=theta<pi 
// costheta=-costheta; 
// -pi<2theta<pi, -pi/2<theta<pi/2 
if ( b<0.0 )// sin2theta<0.0: -pi<2theta<0, -pi/2<theta<0 
sintheta=-sintheta; 
 
// calculate: A=Rt*A*R 
// first save old matrix 
for ( i=0;  i<nn;  i++ ) 
{ 
tip[i]=aa[i*nn+p]; 
tiq[i]=aa[i*nn+q]; 
} 
for ( j=0;  j<nn;  j++ ) 
{ 
tpj[j]=aa[p*nn+j]; 
tqj[j]=aa[q*nn+j]; 
} 
tpp=aa[p*nn+p]; 
tpq=aa[p*nn+q]; 
tqp=aa[q*nn+p]; 
tqq=aa[q*nn+q]; 
// calculate new value 
for ( j=0;  j<nn;  j++ ) 
aa[p*nn+j]=tpj[j]*costheta+tqj[j]*sintheta; 
for ( j=0;  j<nn;  j++ ) 
aa[q*nn+j]=-tpj[j]*sintheta+tqj[j]*costheta; 
for ( i=0;  i<nn;  i++ ) 
aa[i*nn+p]=tip[i]*costheta+tiq[i]*sintheta; 
for ( i=0;  i<nn;  i++ ) 
aa[i*nn+q]=-tip[i]*sintheta+tiq[i]*costheta; 
c=costheta*costheta; 
d=sintheta*sintheta; 
aa[p*nn+p]=tpp*c+tqq*d+tpq*sin2theta; 
aa[q*nn+q]=tpp*d+tqq*c-tpq*sin2theta; 
aa[p*nn+q]=0.0; 
aa[q*nn+p]=0.0; 
 
// calculate S=S*R 
// first save old matrix 
for ( i=0;  i<nn;  i++ ) 
{ 
tip[i]=ss[i*nn+p]; 
tiq[i]=ss[i*nn+q]; 
} 
// calculate new value 
for ( i=0;  i<nn;  i++ ) 
ss[i*nn+p]=tip[i]*costheta+tiq[i]*sintheta; 
for ( i=0;  i<nn;  i++ ) 
ss[i*nn+q]=-tip[i]*sintheta+tiq[i]*costheta; 
} // while ( --count>0 ) 
 
delete []bb; 
} 
 
///////////////////////// 
// methods for CEquationMatrix 
 
BOOL CEquationMatrix::Gauss() 
{ 
double *xx; // solution 
int                 i, j, k, l; 
double a, b, c, d; 
int *ex; // information about column exchange 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( (m_rows+1)!=m_cols ) 
return FALSE; 
 
xx=new double[m_rows]; 
if ( xx==NULL ) 
return FALSE; 
ex=new int[m_rows]; 
if ( ex==NULL ) 
{ 
delete []xx; 
return FALSE; 
} 
 
for ( i=0;  i<m_rows;  i++ ) 
{ 
xx[i]=0.0; 
ex[i]=i; 
} 
 
for ( k=0;  k<m_rows-1;  k++ ) // first m_rows-1 rows 
{ 
d=0.0; 
// search for main element 
for ( i=k;  i<m_rows;  i++ ) 
{ 
for ( j=k;  j<m_cols-1;  j++ ) 
{ 
a=GetElement(i, j); // A[i, j] 
if ( fabs(a)>d )  
{ 
d=fabs(a); 
ex[k]=j; 
l=i; // row 
} 
} 
} 
if ( (d+1.0)==1.0 ) // abnormal 
{ 
delete []ex; 
delete []xx; 
return FALSE; 
} 
 
if ( ex[k]!=k ) // exchange two columns 
{ // A[i, j] <--> A[i, ex[k]] 
for ( i=0;  i<m_rows;  i++ ) 
{ 
a=GetElement(i, ex[k]); 
b=GetElement(i, k); 
SetElement(i, ex[k], b); 
SetElement(i, k, a); 
} 
} 
if ( l!=k ) // exchange two rows 
{ // A[k, j] <--> A[l, j] && B[k] <--> B[l] 
for ( j=k;  j<m_cols;  j++ ) 
{ 
a=GetElement(l, j); 
b=GetElement(k, j); 
SetElement(l, j, b); 
SetElement(k, j, a); 
} 
} 
 
a=GetElement(k, k); // A[k, k] 
SetElement(k, k, 1.0); 
for ( j=k+1;  j<m_cols;  j++ ) 
{ // A[k, j] /= A[k, k] && B[k] /= A[k, k] 
b=GetElement(k, j); 
SetElement(k, j, b/a); 
} 
 
// dissolve 
for ( i=k+1;  i<m_rows;  i++ ) 
{ // A[i, j] -= A[i, k] * A[k, j]/A[k, k]  && A[k, k,]=1 
// B[i] -= A[i, k] * B[k] / A[k, k]  && A[k, k]=1 
a=GetElement(i, k); 
for ( j=k+1;  j<m_cols;  j++ ) 
{ 
b=GetElement(k, j); 
c=GetElement(i, j); 
SetElement(i, j, c-a*b); 
} 
} 
} 
// last row 
d=GetElement(m_rows-1, m_rows-1); 
if ( (d+1.0)==1.0 ) // abnormal 
{ 
delete []ex; 
delete []xx; 
return FALSE; 
} 
a=GetElement(m_rows-1, m_cols-1); 
xx[m_rows-1]=a/d; 
// back iteration 
for ( i=m_rows-2;  i>=0;  i-- ) 
{ 
d=0.0; 
for ( j=i+1;  j<m_rows;  j++ ) 
{ 
b=GetElement(i, j); 
d += b*xx[j]; 
} 
a=GetElement(i, m_cols-1); // B[i] 
xx[i]=a-d; 
} 
// recovery order 
for ( k=m_rows-1;  k>=0;  k-- ) 
{ 
if ( ex[k]!=k ) 
{ 
a=xx[k]; 
xx[k]=xx[ex[k]]; 
xx[ex[k]]=a; 
} 
} 
 
// put resolution to the first column 
for ( i=0;  i<m_rows;  i++ ) 
{ 
SetElement(i, m_cols-1, xx[i]); 
} 
 
delete []ex; 
delete []xx; 
return TRUE; 
} 
 
BOOL CEquationMatrix::Jordan() 
{ 
int                 i, j, k, l; 
double a, b, c, d; 
int *ex; // information about column exchange 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( (m_rows+1)!=m_cols ) 
return FALSE; 
 
ex=new int[m_rows]; 
if ( ex==NULL ) 
return FALSE; 
 
for ( i=0;  i<m_rows;  i++ ) 
ex[i]=i; 
 
for ( k=0;  k<m_rows;  k++ )  
{ 
d=0.0; 
// search for main element 
for ( i=k;  i<m_rows;  i++ ) 
{ 
for ( j=k;  j<m_cols-1;  j++ ) 
{ 
a=GetElement(i, j); // A[i, j] 
if ( fabs(a)>d )  
{ 
d=fabs(a); 
ex[k]=j; 
l=i; // row 
} 
} 
} 
if ( (d+1.0)==1.0 ) // abnormal 
{ 
delete []ex; 
return FALSE; 
} 
 
if ( ex[k]!=k ) // exchange two columns 
{ // A[i, j] <--> A[i, ex[k]] 
for ( i=0;  i<m_rows;  i++ ) 
{ 
a=GetElement(i, ex[k]); 
b=GetElement(i, k); 
SetElement(i, ex[k], b); 
SetElement(i, k, a); 
} 
} 
if ( l!=k ) // exchange two rows 
{ // A[k, j] <--> A[l, j] && B[k] <--> B[l] 
for ( j=k;  j<m_cols;  j++ ) 
{ 
a=GetElement(l, j); 
b=GetElement(k, j); 
SetElement(l, j, b); 
SetElement(k, j, a); 
} 
} 
 
a=GetElement(k, k); // A[k, k] 
SetElement(k, k, 1.0); 
for ( j=k+1;  j<m_cols;  j++ ) 
{ // A[k, j] /= A[k, k] && B[k] /= A[k, k] 
b=GetElement(k, j); 
SetElement(k, j, b/a); 
} 
 
// dissolve 
for ( i=0;  i<m_rows;  i++ ) 
{ // A[i, j] -= A[i, k] * A[k, j]/A[k, k]  && A[k, k,]=1 
// B[i] -= A[i, k] * B[k] / A[k, k]  && A[k, k]=1 
if ( i==k ) 
continue; 
 
a=GetElement(i, k); 
for ( j=k+1;  j<m_cols;  j++ ) 
{ 
b=GetElement(k, j); 
c=GetElement(i, j); 
SetElement(i, j, c-a*b); 
} 
} 
} 
// recovery order 
for ( k=m_rows-1;  k>=0;  k-- ) 
{ 
if ( ex[k]!=k ) 
{ 
a=GetElement(k, m_cols-1); 
b=GetElement(ex[k], m_cols-1); 
SetElement(k, m_cols-1, b); 
SetElement(ex[k], m_cols-1, a); 
} 
} 
 
delete []ex; 
return TRUE; 
} 
 
BOOL CEquationMatrix::Seidel() 
{ 
int                 i, j, l; 
double s, p, t, q; 
double *xx; // solution 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( (m_rows+1)!=m_cols ) 
return FALSE; 
 
xx=new double[m_rows]; 
if ( xx==NULL ) 
return FALSE; 
for ( i=0;  i<m_rows;  i++ ) 
xx[i]=0.0; 
 
// check matrix for diagonal  
for ( i=0;  i<m_rows;  i++ ) 
{ 
s=GetElement(i, i); // A[i, i] 
s=fabs(s); 
p=0.0; 
for ( j=0;  j<m_cols-1;  j++ ) 
{ 
if ( j==i ) 
continue; 
 
t=GetElement(i, j); // A[i, j] 
p += fabs(t); 
} 
if ( p>=s ) 
{ 
delete []xx; 
return FALSE; 
} 
} 
 
l=MATRIXITERATION; 
p=1.0+MATRIXEPSILON; 
while ( p>=MATRIXEPSILON && l>0 ) 
{ 
p=0.0; 
for ( i=0;  i<m_rows;  i++ ) 
{ 
t=xx[i]; 
s=0.0; 
for ( j=0;  j<m_cols-1;  j++ ) 
{ 
if ( j==i ) 
continue; 
s += GetElement(i, j)*xx[j]; 
} 
xx[i]=(GetElement(i, m_rows) - s)/GetElement(i, i); 
 
q=fabs(xx[i]-t); 
if ( q>p ) 
p=q; 
} 
l--; 
} 
for ( i=0;  i<m_rows;  i++ ) 
SetElement(i, m_rows, xx[i]); 
 
delete []xx; 
return TRUE; 
} 
 

 
-- 
※ 来源:．武汉白云黄鹤站WWW bbs.whnet.edu.cn. [FROM: 211.101.187.186]  

华中地区网络中心