求广义特征值的广义雅可比方法
<P><BR>void JacobiEigen(matrix *matA, matrix *matB, double *m_adEigen)<BR>{<BR> // CMatrix matBBak,matR,matRT,matV;<BR> matrix matBBak,matR,matRT,matV;<BR> double dAlfa,dBeta,dBuf,dBuf1;<BR> double dA,dB,dC;<BR> double dError;<BR> int nRow;<BR> int iRow,iCol;<BR> int loop,loop1;<BR> // CArray<double,double&> adBuf;// double向量<BR> // CUIntArray aiBuf; // unsigned int vector<BR> double *adBuf;<BR> unsigned int *aiBuf;</P><P> nRow=matA->m;<BR> MatrixInit(&matBBak,nRow,nRow);<BR> MatrixInit(&matR,nRow,nRow);<BR> MatrixInit(&matRT,nRow,nRow);<BR> MatrixInit(&matV,nRow,nRow);<BR> MatrixCopy(&matBBak,matB);</P>
<P> /* nRow=matA.GetRow();<BR> matBBak.Realloc(nRow,nRow);<BR> matR.Realloc(nRow,nRow);<BR> matRT.Realloc(nRow,nRow);<BR> matV.Realloc(nRow,nRow);<BR> matBBak=matB; */</P>
<P> // matV=0.0;<BR> MatrixZero(&matV);<BR> for(loop=0;loop<nRow;loop++)<BR>matV.a=1.0;<BR>//matV(loop,loop)=1.0;</P>
<P> dError=1.0e-12;</P>
<P> dBuf=0.0;<BR> for(loop=0;loop<nRow;loop++){<BR>for(loop1=loop+1;loop1<nRow;loop1++){<BR> // dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));<BR> dBuf1=p2(matA->a) / (matA->a*matA->a);<BR> if(dBuf<fabs(dBuf1)){<BR> dBuf=fabs(dBuf1);<BR> iRow=loop;iCol=loop1;<BR> }<BR> // dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));<BR> dBuf1=p2(matB->a) / (matB->a*matB->a);<BR> if(dBuf<fabs(dBuf1)){<BR> dBuf=fabs(dBuf1);<BR> iRow=loop;iCol=loop1;<BR> }<BR>}<BR> }</P>
<P> printf("nRow=%ddBuf=%f\n",nRow,dBuf);</P>
<P> while(dBuf>dError){<BR>// matR=0.0;<BR>MatrixZero(&matR);<BR>for(loop=0;loop<nRow;loop++) matR.a=1.0; //matR(loop,loop)=1.0;<BR>// dA=matA(iRow,iRow)*matB(iRow,iCol)-matB(iRow,iRow)*matA(iRow,iCol);<BR>dA=matA->a * matB->a - matB->a * matA->a;<BR>// dC=-matA(iCol,iCol)*matB(iRow,iCol)+matB(iCol,iCol)*matA(iRow,iCol);<BR>dC=-matA->a * matB->a + matB->a * matA->a;<BR>if(dA==0.0&&dC==0.0){<BR> dAlfa=0.0;<BR> // dBeta=-matA(iRow,iCol)/matA(iCol,iCol);<BR> dBeta=-matA->a/matA->a;<BR>}<BR>else{<BR> // dB=matA(iRow,iRow)*matB(iCol,iCol)-matA(iCol,iCol)*matB(iRow,iRow);<BR> dB=matA->a*matB->a - matA->a*matB->a;<BR> if(dB>=0.0){<BR> dBuf=dB/2.0;<BR> dAlfa=-dC/(dBuf+sqrt(dBuf*dBuf-dA*dC));<BR> }<BR> else{<BR> dBuf=dB/2.0;<BR> dAlfa=-dC/(dBuf-sqrt(dBuf*dBuf-dA*dC));<BR> }<BR> dBeta=dA/dC*dAlfa;<BR>}<BR> /* matR(iRow,iCol)=dAlfa;<BR>matR(iCol,iRow)=dBeta;<BR>matRT.Trans(matR);</P>
<P>matA=matRT*matA*matR;<BR>matB=matRT*matB*matR;<BR>matV*=matR;<BR> */<BR>matR.a=dAlfa;<BR>matR.a=dBeta;<BR>Transpose(&matR,&matRT);<BR>// 经过修改后的矩阵相乘功能简单而强大,省去了运算符重载。<BR>MatrixMultiply(&matRT,matA,matA); MatrixMultiply(matA,&matR,matA);<BR>MatrixMultiply(&matRT,matB,matB); MatrixMultiply(matB,&matR,matB);<BR>MatrixMultiply(&matV,&matR,&matV);</P>
<P>printf("A:\n");MatrixPrint(matA); printf("B:\n");MatrixPrint(matB);</P>
<P>dBuf=0.0;<BR>for(loop=0;loop<nRow;loop++){<BR> for(loop1=loop+1;loop1<nRow;loop1++){<BR> // dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));<BR> dBuf1=p2(matA->a)/ (matA->a*matA->a);<BR> if(dBuf<fabs(dBuf1)){<BR> dBuf=fabs(dBuf1);<BR> iRow=loop;iCol=loop1;<BR> }<BR> // dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));<BR> dBuf1=p2(matB->a) / (matB->a*matB->a);<BR> if(dBuf<fabs(dBuf1)){<BR> dBuf=fabs(dBuf1);<BR> iRow=loop;iCol=loop1;<BR> }<BR> }<BR>}<BR>getch();<BR> }</P>
<P> /* matRT.Trans(matV);<BR> matR=matRT*matBBak; */<BR> MatrixTranspose(&matV);<BR> MatrixMultiply(&matRT,&matBBak,&matR);<BR> for(loop=0;loop<nRow;loop++){<BR>dBuf=0.0;<BR>for(loop1=0;loop1<nRow;loop1++){<BR> // dBuf+=matR(loop,loop1)*matV(loop1,loop);<BR> dBuf += matR.a*matV.a;<BR>}<BR>if(dBuf<0.0) dBuf=-dBuf;<BR>dBuf=sqrt(dBuf);<BR>for(loop1=0;loop1<nRow;loop1++)<BR> matV.a /=dBuf;<BR> // matV(loop1,loop)/=dBuf;<BR> }</P>
<P> // adBuf.SetSize(nRow);<BR> adBuf=(double*)malloc(sizeof(double)*nRow);<BR> for(loop=0;loop<nRow;loop++){<BR>// adBuf=matA(loop,loop)/matB(loop,loop);<BR>adBuf=matA->a/matB->a;<BR> }<BR> // aiBuf.SetSize(nRow);<BR> aiBuf=(unsigned int *)malloc(sizeof(unsigned int)*nRow);<BR> aiBuf=0;<BR> for(loop=0;loop<nRow;loop++){<BR>for(loop1=0;loop1<loop;loop1++){<BR> if(fabs(adBuf])>fabs(adBuf)){<BR> // aiBuf.InsertAt(loop1,loop);<BR> for(iRow=nRow-1; iRow>loop1; iRow--)<BR> aiBuf=aiBuf;<BR> aiBuf=loop;<BR> break;<BR> }<BR>}<BR>if(loop1==loop){<BR> aiBuf=loop;<BR>}<BR> }</P>
<P> for(loop=0;loop<nRow;loop++){<BR>m_adEigen=adBuf];<BR> }<BR> for(loop=0;loop<nRow;loop++){<BR>for(loop1=0;loop1<nRow;loop1++){<BR> // matB(loop1,loop)=matV(loop1,aiBuf);<BR> matB->a=matV.a];<BR>}<BR> }</P>
<P> // matrix and vector free<BR> MatrixFree(&matBBak);<BR> MatrixFree(&matR);<BR> MatrixFree(&matRT);<BR> MatrixFree(&matV);<BR> free(adBuf);<BR> free(aiBuf);<BR>}<BR></P>
还有我自编的简单矩阵库
<P>#include <stdio.h><BR>#include <malloc.h><BR>#include <stdlib.h></P><P>#ifndef MATRIXLIB<BR>#define MATRIXLIB</P>
<P>typedef double type;<BR>typedef struct<BR>{<BR> type** a;<BR> int m,n; // m,n for row,col<BR>}matrix;</P>
<P>type** Make2DArray(int row, int col);<BR>void Diliver2DArray(type** a,int row);</P>
<P>type** Make2DArray(int row, int col)<BR>{<BR> type** a;<BR> int i;</P>
<P> if( (a=(type**)malloc(row*sizeof(type*))) == NULL )<BR> {<BR>printf("memory not sufficent!");<BR>exit(0);<BR> }<BR> for( i=0; i<row; i++)<BR>if( (a=(type*)malloc(col*sizeof(type))) == NULL )<BR>{<BR> printf("momory not sifficent!");<BR> exit(0);<BR>}<BR> return a;<BR>}</P>
<P>void Diliver2DArray(type** a,int row)<BR>{<BR> int i;<BR> for( i=0; i<row; i++)<BR>free(a);<BR> free(a);<BR>}</P>
<P>void MatrixInit(matrix* a,int m,int n)<BR>{<BR> a->a = Make2DArray(m,n);<BR> a->m=m; a->n=n;<BR> // 是否应该将元素初始化为0???<BR>}<BR>void MatrixFree(matrix* a)<BR>{<BR> Diliver2DArray(a->a,a->m);<BR>}</P>
<P>// A=0, 矩阵A赋值为0<BR>void MatrixZero(matrix* A)<BR>{<BR> int i,j;<BR> for(i=0; i<A->m; i++)<BR>for(j=0; j<A->n; j++)<BR> A->a=0;<BR>}</P>
<P>// A=B, <BR>void MatrixCopy(matrix *A, matrix *B)<BR>{<BR> int i,j;<BR> if( (A->m!=B->m) || (A->n!=B->n) ) {<BR>printf("cannot add matrix!\n");<BR>return ;<BR> }<BR> for(i=0; i<B->m; i++)<BR>for(j=0; j<B->n; j++)<BR> A->a = B->a;<BR>}</P>
<P>// B=A-1, 矩阵A的逆阵<BR>void Inverse(matrix* A,matrix* B)<BR>{<BR>}</P>
<P>// C=A+B<BR>void MatrixAdd(matrix* A, matrix* B, matrix* C)<BR>{<BR> int i,j;<BR> if( (A->m!=B->m) || (B->m!=C->m) || (A->n!=B->n) || (B->n!=C->n) )<BR> {<BR>printf("cannot add matrix!\n");<BR>return ;<BR> }<BR> for(i=0; i<A->m; i++)<BR>for(j=0; j<A->n; j++)<BR> C->a=A->a+B->a;<BR>}</P>
<P>// c=a*b, 将这个函数功能加强,即 a或b可以跟c相同。<BR>void MatrixMultiply(matrix* a, matrix* b, matrix* c)<BR>{<BR> int m=a->m, n1=a->n, m1=b->m, n=b->n,i,j,k;<BR> double tp=0;<BR> matrix temp;<BR> MatrixInit(&temp,c->m,c->n);<BR> if(m1!=n1)<BR> {<BR>printf("cann't multiply!\n"); exit(0);<BR> }<BR> for(i=0; i<m; i++)<BR>for(j=0; j<n; j++)<BR>{<BR> tp=0;<BR> for(k=0; k<n1; k++)<BR> tp += a->a * b->a;<BR> temp.a=tp;<BR>}<BR> // c=temp;<BR> MatrixCopy(c,&temp);<BR> MatrixFree(&temp);<BR>}</P>
<P>// 方阵A转置<BR>void MatrixTranspose(matrix *A)<BR>{<BR> int i,j;<BR> type temp;<BR> if(A->m!=A->n) {<BR>printf("cannot transpose !");<BR>return ;<BR> }<BR> for(i=0; i<A->m; i++)<BR>for(j=0; j<i; j++) {<BR> temp=A->a;<BR> A->a=A->a;<BR> A->a=temp;<BR>}<BR>}</P>
<P>// B=A的转置<BR>void Transpose(matrix* A, matrix* B)<BR>{<BR> int i,j;<BR> if( (A->m != B->n) || (A->n != B->m) )<BR> {<BR>printf("cannot transpose!\n");<BR>return ;<BR> }<BR> for(i=0; i<B->m; i++)<BR>for(j=0; j<B->n; j++)<BR> B->a=A->a;<BR>}<BR>// A=A*f 矩阵数乘<BR>void MatrixNumber(matrix* A, double f)<BR>{<BR> int i,j;<BR> for(i=0; i<A->m; i++)<BR>for(j=0; j<A->n; j++)<BR> A->a *= f;<BR>}</P>
<P>// c=a*v 矩阵乘向量<BR>void MatrixMulVector(matrix* a, type* v, type* c)<BR>{<BR> int i,j, m=a->m, n=a->n;<BR> double tp=0;<BR> for(i=0; i<m; i++)<BR> {<BR>tp=0;<BR>for(j=0; j<n; j++)<BR> tp += a->a*v;<BR>c=tp;<BR> }<BR>}</P>
<P>// 输出矩阵A,<BR>void MatrixPrint(matrix* A)<BR>{<BR> int i,j;<BR> for(i=0; i<A->m; i++){<BR>for(j=0; j<A->n; j++)<BR> printf("%7.4f ",A->a);<BR>printf("\n");<BR> }<BR>}</P>
<P>#endif</P>
回复:(markailee1)求广义特征值的广义雅可比方法
谢谢分享,希望以后多多支持
页:
[1]