算法及编程语言 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

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

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

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=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) / (matA->a*matA->a); 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) / (matB->a*matB->a); if(dBuf<fabs(dBuf1)){ dBuf=fabs(dBuf1); iRow=loop;iCol=loop1; } } }
 printf("nRow=%ddBuf=%f\n",nRow,dBuf);
 while(dBuf>dError){ // matR=0.0; MatrixZero(&matR); for(loop=0;loop<nRow;loop++) matR.a=1.0; //matR(loop,loop)=1.0; // dA=matA(iRow,iRow)*matB(iRow,iCol)-matB(iRow,iRow)*matA(iRow,iCol); dA=matA->a * matB->a - matB->a * matA->a; // dC=-matA(iCol,iCol)*matB(iRow,iCol)+matB(iCol,iCol)*matA(iRow,iCol); dC=-matA->a * matB->a + matB->a * matA->a; if(dA==0.0&&dC==0.0){ dAlfa=0.0; // dBeta=-matA(iRow,iCol)/matA(iCol,iCol); dBeta=-matA->a/matA->a; } else{ // dB=matA(iRow,iRow)*matB(iCol,iCol)-matA(iCol,iCol)*matB(iRow,iRow); dB=matA->a*matB->a - matA->a*matB->a; 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=dAlfa; matR.a=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)/ (matA->a*matA->a); 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) / (matB->a*matB->a); 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*matV.a; } if(dBuf<0.0) dBuf=-dBuf; dBuf=sqrt(dBuf); for(loop1=0;loop1<nRow;loop1++) matV.a /=dBuf; // matV(loop1,loop)/=dBuf; }
 // adBuf.SetSize(nRow); adBuf=(double*)malloc(sizeof(double)*nRow); for(loop=0;loop<nRow;loop++){ // adBuf=matA(loop,loop)/matB(loop,loop); adBuf=matA->a/matB->a; } // aiBuf.SetSize(nRow); aiBuf=(unsigned int *)malloc(sizeof(unsigned int)*nRow); aiBuf=0; for(loop=0;loop<nRow;loop++){ for(loop1=0;loop1<loop;loop1++){ if(fabs(adBuf])>fabs(adBuf)){ // aiBuf.InsertAt(loop1,loop); for(iRow=nRow-1; iRow>loop1; iRow--) aiBuf=aiBuf; aiBuf=loop; break; } } if(loop1==loop){ aiBuf=loop; } }
 for(loop=0;loop<nRow;loop++){ m_adEigen=adBuf]; } for(loop=0;loop<nRow;loop++){ for(loop1=0;loop1<nRow;loop1++){ // matB(loop1,loop)=matV(loop1,aiBuf); matB->a=matV.a]; } }
 // 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=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 = B->a; }
// 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=A->a+B->a; }
// 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 * b->a; temp.a=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; A->a=A->a; A->a=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=A->a; } // 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 *= 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*v; 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); printf("\n"); } }

#endif

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

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

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

页: [1]

声振论坛's Archiver

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

还有我自编的简单矩阵库

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