求广义特征值的广义雅可比方法

markailee1 · 发表于 2006-6-22 23:30

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

void JacobiEigen(matrix *matA, matrix *matB, double *m_adEigen) { // CMatrix matBBak,matR,matRT,matV; matrix matBBak,matR,matRT,matV; double dAlfa,dBeta,dBuf,dBuf1; double dA,dB,dC; double dError; int nRow; int iRow,iCol; int loop,loop1; // CArray<double,double&> adBuf; // double向量 // CUIntArray aiBuf; // unsigned int vector double *adBuf; unsigned int *aiBuf;
 nRow=matA->m; MatrixInit(&matBBak,nRow,nRow); MatrixInit(&matR,nRow,nRow); MatrixInit(&matRT,nRow,nRow); MatrixInit(&matV,nRow,nRow); MatrixCopy(&matBBak,matB);
 /* nRow=matA.GetRow(); matBBak.Realloc(nRow,nRow); matR.Realloc(nRow,nRow); matRT.Realloc(nRow,nRow); matV.Realloc(nRow,nRow); matBBak=matB; */
 // matV=0.0; MatrixZero(&matV); for(loop=0;loop<nRow;loop++) matV.a[loop][loop]=1.0; //matV(loop,loop)=1.0;
 dError=1.0e-12;
 dBuf=0.0; for(loop=0;loop<nRow;loop++){ for(loop1=loop+1;loop1<nRow;loop1++){ // dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1)); dBuf1=p2(matA->a[loop][loop1]) / (matA->a[loop][loop]*matA->a[loop1][loop1]); if(dBuf<fabs(dBuf1)){ dBuf=fabs(dBuf1); iRow=loop;iCol=loop1; } // dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1)); dBuf1=p2(matB->a[loop][loop1]) / (matB->a[loop][loop]*matB->a[loop1][loop1]); if(dBuf<fabs(dBuf1)){ dBuf=fabs(dBuf1); iRow=loop;iCol=loop1; } } }
 printf("nRow=%d dBuf=%f\n",nRow,dBuf);
 while(dBuf>dError){ // matR=0.0; MatrixZero(&matR); for(loop=0;loop<nRow;loop++) matR.a[loop][loop]=1.0; //matR(loop,loop)=1.0; // dA=matA(iRow,iRow)*matB(iRow,iCol)-matB(iRow,iRow)*matA(iRow,iCol); dA=matA->a[iRow][iRow] * matB->a[iRow][iCol] - matB->a[iRow][iRow] * matA->a[iRow][iCol]; // dC=-matA(iCol,iCol)*matB(iRow,iCol)+matB(iCol,iCol)*matA(iRow,iCol); dC=-matA->a[iCol][iCol] * matB->a[iRow][iCol] + matB->a[iCol][iCol] * matA->a[iRow][iCol]; if(dA==0.0&&dC==0.0){ dAlfa=0.0; // dBeta=-matA(iRow,iCol)/matA(iCol,iCol); dBeta=-matA->a[iRow][iCol]/matA->a[iCol][iCol]; } else{ // dB=matA(iRow,iRow)*matB(iCol,iCol)-matA(iCol,iCol)*matB(iRow,iRow); dB=matA->a[iRow][iRow]*matB->a[iCol][iCol] - matA->a[iCol][iCol]*matB->a[iRow][iRow]; if(dB>=0.0){ dBuf=dB/2.0; dAlfa=-dC/(dBuf+sqrt(dBuf*dBuf-dA*dC)); } else{ dBuf=dB/2.0; dAlfa=-dC/(dBuf-sqrt(dBuf*dBuf-dA*dC)); } dBeta=dA/dC*dAlfa; } /* matR(iRow,iCol)=dAlfa; matR(iCol,iRow)=dBeta; matRT.Trans(matR);
 matA=matRT*matA*matR; matB=matRT*matB*matR; matV*=matR; */ matR.a[iRow][iCol]=dAlfa; matR.a[iCol][iRow]=dBeta; Transpose(&matR,&matRT); // 经过修改后的矩阵相乘功能简单而强大，省去了运算符重载。 MatrixMultiply(&matRT,matA,matA); MatrixMultiply(matA,&matR,matA); MatrixMultiply(&matRT,matB,matB); MatrixMultiply(matB,&matR,matB); MatrixMultiply(&matV,&matR,&matV);
 printf("A:\n");MatrixPrint(matA); printf("B:\n");MatrixPrint(matB);
 dBuf=0.0; for(loop=0;loop<nRow;loop++){ for(loop1=loop+1;loop1<nRow;loop1++){ // dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1)); dBuf1=p2(matA->a[loop][loop1])/ (matA->a[loop][loop]*matA->a[loop1][loop1]); if(dBuf<fabs(dBuf1)){ dBuf=fabs(dBuf1); iRow=loop;iCol=loop1; } // dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1)); dBuf1=p2(matB->a[loop][loop1]) / (matB->a[loop][loop]*matB->a[loop1][loop1]); if(dBuf<fabs(dBuf1)){ dBuf=fabs(dBuf1); iRow=loop;iCol=loop1; } } } getch(); }
 /* matRT.Trans(matV); matR=matRT*matBBak; */ MatrixTranspose(&matV); MatrixMultiply(&matRT,&matBBak,&matR); for(loop=0;loop<nRow;loop++){ dBuf=0.0; for(loop1=0;loop1<nRow;loop1++){ // dBuf+=matR(loop,loop1)*matV(loop1,loop); dBuf += matR.a[loop][loop1]*matV.a[loop1][loop]; } if(dBuf<0.0) dBuf=-dBuf; dBuf=sqrt(dBuf); for(loop1=0;loop1<nRow;loop1++) matV.a[loop1][loop] /=dBuf; // matV(loop1,loop)/=dBuf; }
 // adBuf.SetSize(nRow); adBuf=(double*)malloc(sizeof(double)*nRow); for(loop=0;loop<nRow;loop++){ // adBuf[loop]=matA(loop,loop)/matB(loop,loop); adBuf[loop]=matA->a[loop][loop]/matB->a[loop][loop]; } // aiBuf.SetSize(nRow); aiBuf=(unsigned int *)malloc(sizeof(unsigned int)*nRow); aiBuf[0]=0; for(loop=0;loop<nRow;loop++){ for(loop1=0;loop1<loop;loop1++){ if(fabs(adBuf[aiBuf[loop1]])>fabs(adBuf[loop])){ // aiBuf.InsertAt(loop1,loop); for(iRow=nRow-1; iRow>loop1; iRow--) aiBuf[iRow]=aiBuf[iRow-1]; aiBuf[loop1]=loop; break; } } if(loop1==loop){ aiBuf[loop1]=loop; } }
 for(loop=0;loop<nRow;loop++){ m_adEigen[loop]=adBuf[aiBuf[loop]]; } for(loop=0;loop<nRow;loop++){ for(loop1=0;loop1<nRow;loop1++){ // matB(loop1,loop)=matV(loop1,aiBuf[loop]); matB->a[loop1][loop]=matV.a[loop1][aiBuf[loop]]; } }
 // matrix and vector free MatrixFree(&matBBak); MatrixFree(&matR); MatrixFree(&matRT); MatrixFree(&matV); free(adBuf); free(aiBuf); }

markailee1 · 发表于 2006-6-22 23:32

#include <stdio.h> #include <malloc.h> #include <stdlib.h>
#ifndef MATRIXLIB #define MATRIXLIB
typedef double type; typedef struct { type** a; int m,n; // m,n for row,col }matrix;
type** Make2DArray(int row, int col); void Diliver2DArray(type** a,int row);
type** Make2DArray(int row, int col) { type** a; int i;
 if( (a=(type**)malloc(row*sizeof(type*))) == NULL ) { printf("memory not sufficent!"); exit(0); } for( i=0; i<row; i++) if( (a=(type*)malloc(col*sizeof(type))) == NULL ) { printf("momory not sifficent!"); exit(0); } return a; }
void Diliver2DArray(type** a,int row) { int i; for( i=0; i<row; i++) free(a); free(a); }
void MatrixInit(matrix* a,int m,int n) { a->a = Make2DArray(m,n); a->m=m; a->n=n; // 是否应该将元素初始化为0??? } void MatrixFree(matrix* a) { Diliver2DArray(a->a,a->m); }
// A=0, 矩阵A赋值为0 void MatrixZero(matrix* A) { int i,j; for(i=0; i<A->m; i++) for(j=0; j<A->n; j++) A->a[j]=0; }
// A=B, void MatrixCopy(matrix *A, matrix *B) { int i,j; if( (A->m!=B->m) || (A->n!=B->n) ) { printf("cannot add matrix!\n"); return ; } for(i=0; i<B->m; i++) for(j=0; j<B->n; j++) A->a[j] = B->a[j]; }
// B=A-1, 矩阵A的逆阵 void Inverse(matrix* A,matrix* B) { }
// C=A+B void MatrixAdd(matrix* A, matrix* B, matrix* C) { int i,j; if( (A->m!=B->m) || (B->m!=C->m) || (A->n!=B->n) || (B->n!=C->n) ) { printf("cannot add matrix!\n"); return ; } for(i=0; i<A->m; i++) for(j=0; j<A->n; j++) C->a[j]=A->a[j]+B->a[j]; }
// c=a*b，将这个函数功能加强，即 a或b可以跟c相同。 void MatrixMultiply(matrix* a, matrix* b, matrix* c) { int m=a->m, n1=a->n, m1=b->m, n=b->n, i,j,k; double tp=0; matrix temp; MatrixInit(&temp,c->m,c->n); if(m1!=n1) { printf("cann't multiply!\n"); exit(0); } for(i=0; i<m; i++) for(j=0; j<n; j++) { tp=0; for(k=0; k<n1; k++) tp += a->a[k] * b->a[k][j]; temp.a[j]=tp; } // c=temp; MatrixCopy(c,&temp); MatrixFree(&temp); }
// 方阵A转置 void MatrixTranspose(matrix *A) { int i,j; type temp; if(A->m!=A->n) { printf("cannot transpose !"); return ; } for(i=0; i<A->m; i++) for(j=0; j<i; j++) { temp=A->a[j]; A->a[j]=A->a[j]; A->a[j]=temp; } }

// B=A的转置 void Transpose(matrix* A, matrix* B) { int i,j; if( (A->m != B->n) || (A->n != B->m) ) { printf("cannot transpose!\n"); return ; } for(i=0; i<B->m; i++) for(j=0; j<B->n; j++) B->a[j]=A->a[j]; } // A=A*f 矩阵数乘 void MatrixNumber(matrix* A, double f) { int i,j; for(i=0; i<A->m; i++) for(j=0; j<A->n; j++) A->a[j] *= f; }
// c=a*v 矩阵乘向量 void MatrixMulVector(matrix* a, type* v, type* c) { int i,j, m=a->m, n=a->n; double tp=0; for(i=0; i<m; i++) { tp=0; for(j=0; j<n; j++) tp += a->a[j]*v[j]; c=tp; } }
// 输出矩阵A, void MatrixPrint(matrix* A) { int i,j; for(i=0; i<A->m; i++){ for(j=0; j<A->n; j++) printf("%7.4f ",A->a[j]); printf("\n"); } }

#endif

风花雪月 · 发表于 2006-6-25 08:57

谢谢分享，希望以后多多支持

账号		自动登录	找回密码
密码			我要加入

[经典算法] 求广义特征值的广义雅可比方法

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

还有我自编的简单矩阵库

回复：（markailee1）求广义特征值的广义雅可比方法