最初に

定義式

プログラムリスト

#define C 3			// column
#define R 3			// rank
double factors[C * R]; //行列の定義			//この行はで宣言
void MultipleAbove(CMatrix* a, int seigen = 0);	//この行はで宣言
void MultipleBehind(CMatrix* a, int seigen = 0);	//この行はで宣言
BOOL HouseHolder(double *NORM, double *det);	//この行はで宣言

//ハウスホルダー変換を利用した直交化
BOOL CMatrix::HouseHolder(double *NORM, double *det)
{
	register double len, sigma;
	register int i, j, k, l, cr;
	cr = C * R;
	CMatrix* Q= new CMatrix(R, R);
	CMatrix* A= new CMatrix(C, R, factors);
	double* t= new double[R];

	for(k = 0; k < R; k++){

		//(a, 0, 0, ･･･, 0)化するベクトル t を計算
		for(l = R-k, len = 0.0, i=0; i < l; i++)
			t[i] = A->factors[(k + i) * R + k], len += t[i] * t[i];

		//lenとt[0]は同符号にして桁落ちを回避する
		len = (t[0] < 0)? -sqrt(len): sqrt(len);
		sigma = len, t[0] += len,
		len *= t[0];			//len = sqrt(len * (len + A->factors[k * R + k]) * 2.0);
						//sqrtを使わないであとで手直しをする

		if(len == 0.0){
			if(k == R - 1) break;	//最後の行ではlen==0は操作の必要なしを表す
			delete Q;
			delete A;
			delete [] t;
			return FALSE;
		}

		//Q'を計算
		//Q->unitmatrix(R);
		for(i = 0; i < R; i++){
			for(j = 0; j < R; j++){
				Q->factors[i * R + j] = 0.0;
			}
			Q->factors[i * (R + 1)] = 1.0;
		}

		for(i = 0; i < l; i++){
			for(j = 0; j < l; j++){
				Q->factors[(k + i) * R + k + j] -= t[i] * t[j] /len; //lenで割ってt[i],t[j]を単位行列化する
			}
		}

		//ΠQ×Aの計算(理論的に  [a11  z~])
		//ΠQ×Aの計算(理論的に  [ y   B ])

		//AにQを前から掛ける
		A->MultipleAbove(Q);

		//一応零に近くなるが，0.0を代入して保証する
		for(i = k + 1; i < C; i++){
			A->factors[i * R + k] = 0.0;
		}

		//ΠQの計算
		if(k>0){
			//後ろからQをかける
			MultipleBehind(Q);
		}else{
			for(i = 0;i < cr; i++){
				factors[i] = Q->factors[i];
			}
		}
	}

	//直交行列のノルムの計算
	if(NORM != NULL){
		for(k = 0; k < R; k++){
			NORM[k] = A->factors[k * (R + 1)] * A->factors[k * (R + 1)];
		}
	}

	//元の行列のデターミナントの計算
	if(det ! =NULL){
		for(*det = 1.0, k = 0; k < R; k++){
			*det *= A->factors[k * (R + 1)];
		}
	}

	delete Q;
	delete A;
	delete [] t;

	return TRUE;
}

//(*this) = (a) * (*this)を計算する //seigen:行数制限（いらないデータは計算しない）
void CMatrix::MultipleAbove(CMatrix* a, int seigen) //前からaを掛ける
{
	if(a->R != C){
		AfxMessageBox("CMatrix operator multiple error");
		exit(0);
	}

	int i, j, k;
	double sam;
	double* stack;
	stack = new double [C * R]; //行と列を入れ替えたスタックを用意する

	for(i = 0; i < C; i++){
		for(j = 0; j < R; j++){
			stack[j * C + i]=factors[i * R + j];
		}
	}

	if(C != a-> C){
		//新しい大きさのマトリクスに変える
		delete [] factors;
		C = a->C;
		factors = new double[C * R];
	}

	if(seigen == 0) seigen = C;			//seigen==0は制限なしの意味

	for(i = 0; i < seigen; i++){
		for(j = 0; j < R; j++){
			for(sam = 0.0, k = 0; k < C; k++){
				sam += a->factors[i * a-> R + k]*stack[j * C + k];
			}
			factors[i * R + j] = sam;
		}
	}
	delete [] stack;
}

//(*this)=(*this)*(a)を計算する //seigen:行数制限（いらないデータは計算しない）
void CMatrix::MultipleBehind(CMatrix* a, int seigen) //後ろからaを掛ける
{
	if(R != a->C){
		AfxMessageBox("CMatrix operator multiple error");
		exit(0);
	}

	int i, j, k;
	double sam;
	double* stack;
	stack = new double [C * R];

	for(i = 0;i < C * R; i++){
		stack[i] = factors[i];
	}

	if(R != a->R){
		//新しい大きさのマトリクスに変える
		delete [] factors;
		R = a->R;
		factors = new double[C * R];
	}

	if(seigen == 0) seigen = C;			//seigen==0は制限なしの意味

	for(i = 0; i < seigen; i++){
		for(j = 0;j < R; j++){
			for(sam = 0.0, k = 0; k < R; k++){
				sam += stack[i * R + k] * a-> factors[k * a->R + j];
			}
			factors[i * R + j] = sam;
		}
	}
	delete [] stack;
}