Bidiagonalization using Householder transformations

The previous post was a discussion on employing Householder transformations to perform a QR decomposition. This post will be short. I've had this code lying around for a while now and thought I would make it available. The process of bidiagonalization using Householder transformations amounts to nothing more than alternating by left and right transformations.

The cMatrix::householderBidiagonalization() method:

void cMatrix::householderBidiagonalization(cMatrix& Q, cMatrix& R, cMatrix& S) {
	double mag, alpha;

	cMatrix u(m, 1),  v(m, 1),
		u_(n, 1), v_(n, 1);

	cMatrix P(m, m),  I(m, m),
		P_(n, n), I_(n, n);

	Q = cMatrix(m, m);
	R = *this;
	S = cMatrix(n, n);

	for (int i = 0; i < n; i++) {
		u.zero(); v.zero();
		
		mag = 0.0;
		for (int j = i; j < m; j++) {
			u.A[j] = R.A[j * n + i];
			mag += u.A[j] * u.A[j];
		}
		mag = sqrt(mag);
		
		alpha = u.A[i] < 0 ? mag : -mag;

		mag = 0.0;
		for (int j = i; j < m; j++) {
			v.A[j] = j == i ? u.A[j] + alpha : u.A[j];
			mag += v.A[j] * v.A[j];
		}
		mag = sqrt(mag);

		if (mag > 0.0000000001) {
			for (int j = i; j < m; j++) v.A[j] /= mag;

			P = I - (v * v.transpose()) * 2.0;

			R = P * R;
			Q = Q * P;
		}
/////////////////////////
		u_.zero(); v_.zero();

		mag = 0.0;
		for (int j = i + 1; j < n; j++) {
			u_.A[j] = R.A[i * n + j];
			mag += u_.A[j] * u_.A[j];
		}

		mag = sqrt(mag);

		alpha = u_.A[i + 1] < 0 ? mag : -mag;

		mag = 0.0;
		for (int j = i + 1; j < n; j++) {
			v_.A[j] = j == i + 1 ? u_.A[j] + alpha : u_.A[j];
			mag += v_.A[j] * v_.A[j];
		}
		mag = sqrt(mag);

		if (mag > 0.0000000001) {
			for (int j = i + 1; j < n; j++) v_.A[j] /= mag;

			P_ = I_ - (v_ * v_.transpose()) * 2.0;

			R = R * P_;
			S = P_ * S;
		}
	}
}

Download the source: qr_householder_bidiagonalization.cc.bz2