matrix3.txt

几个算法程序

武汉白云黄鹤站∶精华区武汉白云黄鹤站∶精华区
发信人: sweeping (simon), 信区: Algorithm WWW-POST 
标  题: 一些矩阵源码3 
发信站: 武汉白云黄鹤站 (Thu Nov  8 17:42:10 2001) , 转信 
 
//////////////////////////// 
// methods for CSquareMatrix 
 
void CSquareMatrix::SetToIdentity() 
{ 
int i, n; 
 
if ( m_rows<1 || m_cols<1 ) 
return; 
if ( m_rows!=m_cols ) // not an nxn matrix 
return; 
 
n=m_rows*m_cols; 
for ( i=0;  i<n;  i++ ) 
m_a[i]=0.0; 
for ( i=0;  i<m_rows;  i++ ) 
m_a[i*m_cols+i]=1.0; 
} 
 
double CSquareMatrix::Determinantal() 
{ 
int                 i, j, k, l; 
double a, b, c, d; 
int *ex; // information about column exchange 
int sign; 
 
if ( m_rows<1 || m_cols<1 ) 
return 0.0; 
if ( m_rows!=m_cols ) 
return 0.0; 
 
ex=new int[m_rows]; 
if ( ex==NULL ) 
return 0.0; 
 
for ( i=0;  i<m_rows;  i++ ) 
ex[i]=i; 
sign=1; 
 
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; 
return 0.0; 
} 
 
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); 
} 
sign=-sign; 
} 
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); 
} 
sign=-sign; 
} 
 
a=GetElement(k, k); // A[k, k] 
// dissolve 
for ( i=k+1;  i<m_rows;  i++ ) 
{ // A[i, j] -= A[i, k] * A[k, j]/A[k, k]   
d=GetElement(i, k); 
for ( j=k+1;  j<m_cols;  j++ ) 
{ 
b=GetElement(k, j); 
c=GetElement(i, j); 
SetElement(i, j, c-d*b/a); 
} 
} 
} 
// last row 
d=GetElement(m_rows-1, m_rows-1); 
if ( (d+1.0)==1.0 ) // abnormal 
{ 
delete []ex; 
return 0.0; 
} 
 
// put resolution to the first column 
a=1.0; 
for ( i=0;  i<m_rows;  i++ ) 
a *= GetElement(i, i); 
if ( sign<0 ) 
a = -a; 
 
delete []ex; 
return a; 
} 
 
BOOL CSquareMatrix::Inverse() 
{ 
int                 i, j, k; 
double a, b, c, d; 
int *colex; // information about column exchange 
int *rowex; // information about row exchange 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( m_rows!=m_cols ) 
return FALSE; 
 
colex=new int[m_cols]; 
if ( colex==NULL ) 
return FALSE; 
rowex=new int[m_rows]; 
if ( rowex==NULL ) 
{ 
delete []colex; 
return FALSE; 
} 
 
for ( i=0;  i<m_rows;  i++ ) 
{ 
rowex[i]=i; 
colex[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;  j++ ) 
{ 
a=GetElement(i, j); // A[i, j] 
if ( fabs(a)>d )  
{ 
d=fabs(a); 
rowex[k]=i; 
colex[k]=j; 
} 
} 
} 
if ( (d+1.0)==1.0 ) // abnormal 
{ 
delete []rowex; 
delete []colex; 
return FALSE; 
} 
 
if ( rowex[k]!=k ) // exchange two rows 
{ // A[k, j] <--> A[rowex[k], j] 
for ( j=0;  j<m_cols;  j++ ) 
{ 
a=GetElement(rowex[k], j); 
b=GetElement(k, j); 
SetElement(rowex[k], j, b); 
SetElement(k, j, a); 
} 
} 
if ( colex[k]!=k ) // exchange two columns 
{ // A[i, j] <--> A[i, colex[k]] 
for ( i=0;  i<m_rows;  i++ ) 
{ 
a=GetElement(i, colex[k]); 
b=GetElement(i, k); 
SetElement(i, colex[k], b); 
SetElement(i, k, a); 
} 
} 
 
a=GetElement(k, k); // A[k, k]=1/A[k, k] 
SetElement(k, k, 1.0/a); 
 
a=GetElement(k, k);  
for ( j=0;  j<m_cols;  j++ ) 
{ // A[k, j] = A[k, j] * A[k, k] 
if ( j==k ) 
continue; 
 
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] 
if ( i==k ) 
continue; 
 
d=GetElement(i, k); 
for ( j=0;  j<m_cols;  j++ ) 
{ 
if ( j==k ) 
continue; 
 
b=GetElement(k, j); 
c=GetElement(i, j); 
SetElement(i, j, c-d*b); 
} 
} 
 
a=GetElement(k, k);  
for ( i=0;  i<m_rows;  i++ ) 
{ // A[i, k] = A[i, k] * A[k, k] 
if ( i==k ) 
continue; 
 
b=GetElement(i, k); 
SetElement(i, k, -b*a); 
} 
} 
// recovery order 
for ( k=m_rows-1;  k>=0;  k-- ) 
{ 
for ( j=0;  j<m_cols;  j++ ) 
{ 
if ( colex[k]!=k ) 
{ 
a=GetElement(k, j); 
b=GetElement(colex[k], j); 
SetElement(k, j, b); 
SetElement(colex[k], j, a); 
} 
} 
for ( i=0;  i<m_rows;  i++ ) 
{ 
if ( rowex[k]!=k ) 
{ 
a=GetElement(i, k); 
b=GetElement(i, rowex[k]); 
SetElement(i, k, b); 
SetElement(i, rowex[k], a); 
} 
} 
} 
 
delete []rowex; 
delete []colex; 
return TRUE; 
} 
 
BOOL CSquareMatrix::LUDecompose(CSquareMatrix *pL, CSquareMatrix *pU) 
{ 
int                 i, j, k; 
double a, b, c, d; 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( m_rows!=m_cols ) 
return FALSE; 
if ( pL==NULL || pU==NULL ) 
return FALSE; 
if ( pL->GetRows()!=m_rows || pL->GetRows()!=pL->GetColumns() ) 
if ( !pL->Resize(m_rows) ) 
return FALSE; 
if ( pU->GetRows()!=m_rows || pU->GetRows()!=pU->GetColumns() ) 
if ( !pU->Resize(m_rows) ) 
return FALSE; 
 
for ( k=0;  k<m_rows-1;  k++ )  
{ 
d=GetElement(k, k); // A[k, k] 
if ( (d+1.0)==1.0 ) // abnormal 
return FALSE; 
 
for ( i=k+1;  i<m_cols;  i++ )  
{ // A[i, k] /= A[k, k] 
a=GetElement(i, k);  
SetElement(i, k, a/d);  
} 
 
for ( i=k+1;  i<m_rows;  i++ ) 
{ // A[i, j] -= A[i, k] * A[k, j] 
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); 
} 
} 
} 
// decompose 
for ( i=0;  i<m_rows;  i++ ) 
{ 
for ( j=0;  j<i;  j++ ) 
{ 
a=GetElement(i, j); 
pL->SetElement(i, j, a); 
pU->SetElement(i, j, 0.0); 
} 
 
pL->SetElement(i, i, 1.0); 
a=GetElement(i, i); 
pU->SetElement(i, j, a); 
 
for ( j=i+1;  j<m_cols;  j++ ) 
{ 
pL->SetElement(i, j, 0.0); 
a=GetElement(i, j); 
pU->SetElement(i, j, a); 
} 
} 
 
return TRUE; 
} 
 
BOOL CSquareMatrix::QRDecompose(CSquareMatrix *pQ) 
{ 
int                 i, j, k, l; 
double a, b; 
double u, t; 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( m_rows!=m_cols ) 
return FALSE; 
if ( pQ==NULL ) 
return FALSE; 
if ( pQ->GetRows()!=m_rows || pQ->GetRows()!=pQ->GetColumns() ) 
if ( !pQ->Resize(m_rows) ) 
return FALSE; 
 
pQ->SetToIdentity(); 
 
for ( k=0;  k<m_rows-1;  k++ )  
{ 
u=0.0; 
for ( i=k;  i<m_rows;  i++ ) 
{ 
a=GetElement(i, k);  
if ( fabs(a)>u ) 
u=fabs(a); 
} 
 
b=0.0; 
for ( i=k;  i<m_rows;  i++ ) 
{ 
t=GetElement(i, k);  
t /= u; 
b += t*t; 
} 
 
a=GetElement(k, k); 
if ( a>0.0 ) 
u = -u; 
b=u*sqrt(b); 
 
if ( (b+1.0)==1.0 ) // abnormal 
return FALSE; 
 
a=GetElement(k, k); 
u=sqrt(2.0*b*(b-a)); 
if ( (u+1.0)!=1.0 ) // normal 
{ 
a=GetElement(k, k); 
a=(a-b)/u; 
SetElement(k, k, a); // A[k, k]=(A[k, k]-b)/u 
for ( i=k+1;  i<m_rows;  i++ ) 
{ // A[i, k]/=u 
a=GetElement(i, k); 
SetElement(i, k, a/u); 
} 
for ( j=0;  j<m_cols;  j++ ) 
{ 
t=0.0; 
for ( l=k;  l<m_rows;  l++ ) 
t += GetElement(l, k)*pQ->GetElement(l, j); 
for ( i=k;  i<m_rows;  i++ ) 
{ // Q[i, j] -= 2*t*A[i, k] 
a=pQ->GetElement(i, j); 
a -= 2.0*t*GetElement(i, k); 
pQ->SetElement(i, j, a); 
} 
} 
for ( j=k+1;  j<m_cols;  j++ ) 
{ 
t=0.0; 
for ( l=k;  l<m_rows;  l++ ) 
t += GetElement(l, k)*GetElement(l, j); 
for ( i=k;  i<m_rows;  i++ ) 
{ // A[i, j] -= 2*t*A[i, k] 
a=GetElement(i, j); 
a -= 2.0*t*GetElement(i, k); 
SetElement(i, j, a); 
} 
SetElement(k, k, b); // A[k, k]=b 
} 
for ( i=k+1;  i<m_rows;  i++ ) // A[i, k]=0 
SetElement(i, k, 0.0); 
} 
} 
 
// transpose 
for ( i=0;  i<m_rows-1;  i++ ) 
{ // Q[i, j] = Q[j, i] 
for ( j=i+1;  j<m_cols;  j++ ) 
{ 
a=pQ->GetElement(i, j); 
b=pQ->GetElement(j, i); 
pQ->SetElement(i, j, b); 
pQ->SetElement(j, i, a); 
} 
} 
 
return TRUE; 
} 
 
BOOL CSquareMatrix::Hessenberg() 
{ 
int                 i, j, k; 
double a, b, c, d; 
 
if ( m_rows<1 || m_cols<1 ) 
return FALSE; 
if ( m_rows!=m_cols ) 
return FALSE; 
 
for ( k=1;  k<m_rows-1;  k++ )  
{ 
// select main element 
d=0.0; 
i=k; 
for ( j=k;  j<m_cols;  j++ ) 
{ 
a=GetElement(j, k-1);  
if ( fabs(a)>fabs(d) ) 
{ 
d=a; 
i=j; 
} 
} 
 
// exchange 
if ( i!=k ) 
{ 
for ( j=k-1;  j<m_cols;  j++ ) 
{ // A[i, j] <--> A[k, j] 
a=GetElement(i, j);  
b=GetElement(k, j);  
SetElement(k, j, a);  
SetElement(i, j, b);  
} 
for ( j=0;  j<m_rows;  j++ ) 
{ // A[j, i] <--> A[j, k] 
a=GetElement(j, i);  
b=GetElement(j, k);  
SetElement(j, k, a);  
SetElement(j, i, b);  
} 
} 
 
// dissolve 
if ( (d+1.0)!=1.0 ) // normal 
{ 
for ( i=k+1;  i<m_rows;  i++ ) 
{ 
a=GetElement(i, k-1)/d; 
SetElement(i, k-1, 0.0); 
 
for ( j=k;  j<m_cols;  j++ ) 
{ // A[i, j] -= a*A [k, j] 
b=GetElement(i, j); 
c=GetElement(k, j); 
SetElement(i, j, b-c*a); 
} 
 
for ( j=1;  j<m_rows;  j++ ) 
{ // A[j, k] += a*A [j, i] 
b=GetElement(j, k); 
c=GetElement(j, i); 
SetElement(j, k, b+c*a); 
} 
} 
} 
} 
return TRUE; 
} 
 
CMatrix *CSquareMatrix::QREigenvalue() 
{ 
CMatrix *pUV; 
int j, k, l, m, s, t; 
double a, b, c, d, e, f, w, p, q, r, x, y, z, ss; 
double sgn; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
if ( m_rows!=m_cols ) 
return NULL; 
 
pUV=new CMatrix(m_rows, 2); 
if ( pUV==NULL ) 
return NULL; 
pUV->SetAllElements(0.0); 
 
if ( !Hessenberg() ) 
{ 
delete pUV; 
return NULL; 
} 
 
m=m_rows-1; 
t=0; // iteration times 
while ( m!=-1 ) 
{ 
l=m; 
while ( l>0 ) 
{ 
a=GetElement(l, l-1); 
b=GetElement(l-1, l-1); 
c=GetElement(l, l); 
if ( fabs(a)<=(MATRIXEPSILON*(fabs(b)+fabs(c))) ) 
break; 
l--; 
} 
if ( l==m ) // one order block 
{ 
pUV->SetElement(m, 0, GetElement(m, m)); 
pUV->SetElement(m, 1, 0.0); 
m--; 
t=0; 
} 
else 
{ 
if ( l==(m-1) ) // two order block 
{ 
b=-(GetElement(m, m)+GetElement(m-1, m-1)); 
c=GetElement(m, m)*GetElement(m-1, m-1) 
-GetElement(m, m-1)*GetElement(m-1, m); 
w=b*b-4.0*c; 
y=sqrt(fabs(w)); 
if ( w>0.0 ) 
{ 
sgn=0.0; 
if ( b>0.0 ) 
sgn=1.0; 
else if ( b<0.0 ) 
sgn=-1.0; 
a=(-b-sgn*y)/2.0; 
pUV->SetElement(m, 0, a); 
a=pUV->GetElement(m, 0); 
pUV->SetElement(m-1, 0, c/a); 
pUV->SetElement(m, 1, 0.0); 
pUV->SetElement(m-1, 1, 0.0); 
} 
else 
{ 
a=-b/2.0; 
pUV->SetElement(m, 0, a); 
pUV->SetElement(m-1, 0, a); 
pUV->SetElement(m, 1, y/2.0); 
pUV->SetElement(m-1, 1, -y/2.0); 
} 
m-=2; 
t=0; 
} 
else // if l==m-1 
{ 
if ( t>=MATRIXITERATION ) 
{ 
delete pUV; 
return NULL; 
} 
 
t++; 
for ( j=l+2;  j<=m;  j++ ) 
SetElement(j, j-2, 0.0); 
for ( j=l+3;  j<=m;  j++ ) 
SetElement(j, j-3, 0.0); 
for ( k=l;  k<m;  k++ ) 
{ 
if ( k!=l ) 
{ 
p=GetElement(k, k-1); 
q=GetElement(k+1, k-1); 
r=0.0; 
if ( k!=(m-1) ) 
r=GetElement(k+2, k-1); 
} 
else 
{ 
x=GetElement(m, m)+GetElement(m-1, m-1); 
y=GetElement(m-1, m-1)*GetElement(m, m) 
-GetElement(m-1, m)*GetElement(m, m-1); 
p=GetElement(l, l)*(GetElement(l, l)-x) 
+GetElement(l, l+1)*GetElement(l+1, l)+y; 
q=GetElement(l+1, l)*(GetElement(l, l)+GetElement(l+1, l+1)-x); 
r=GetElement(l+1, l)*GetElement(l+2, l+1); 
} 
sgn=0.0; 
if ( p>0.0 ) 
sgn=1.0; 
else if ( p<0.0 ) 
sgn=-1.0; 
ss=sgn*sqrt(p*p+q*q+r*r); 
if ( (fabs(ss)+1.0)!=1.0 ) 
{ 
if ( k!=l ) 
SetElement(k, k-1, -ss); 
d=-q/ss; 
e=-r/ss; 
x=-p/ss; 
y=-x-e*r/(p+ss); 
f=d*r/(p+ss); 
z=-x-d*q/(p+ss); 
for ( j=k;  j<=m;  j++ ) 
{ 
p=x*GetElement(k, j)+d*GetElement(k+1, j); 
q=d*GetElement(k, j)+y*GetElement(k+1, j); 
r=e*GetElement(k, j)+f*GetElement(k+1, j); 
if ( k!=(m-1) ) 
{ 
p +=e*GetElement(k+2, j); 
q +=f*GetElement(k+2, j); 
a = r+z*GetElement(k+2, j); 
SetElement(k+2, j, a); 
} 
SetElement(k+1, j, q); 
SetElement(k, j, p); 
} 
j=k+3; 
if ( j>=m ) 
j=m; 
for ( s=l;  s<=j;  s++ ) 
{ 
p=x*GetElement(s, k)+d*GetElement(s, k+1); 
q=d*GetElement(s, k)+y*GetElement(s, k+1); 
r=e*GetElement(s, k)+f*GetElement(s, k+1); 
if ( k!=(m-1) ) 
{ 
p +=e*GetElement(s, k+2); 
q +=f*GetElement(s, k+2); 
a = r+z*GetElement(s, k+2); 
SetElement(s, k+2, a); 
} 
SetElement(s, k+1, q); 
SetElement(s, k, p); 
} 
} 
} 
} 
} 
} 
return pUV; 
} 

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

华中地区网络中心