Function UdUrecompose_transpose

bayes/src/UdU.cpp:602–644

void UdUrecompose_transpose (RowMatrix& M)
/* In-place recomposition of Symmetric matrix from U'dU factor store in UD format
 *  Generally used for recomposing result of UdUinverse
 * Note definiteness of result depends purely on diagonal(M)
 *  i.e. if d is positive definite (>0) then result is positive definite
 * Reference: A+G p.223
 * In place computation uses simple structure of solution due to triangular zero elements
 *  Defn: R = (U' d) row i , C = U column j   -> M(i,j) = R dot C;
 *  However M(i,j) only dependent R(k<=i), C(k<=j) due to zeros
 *  Therefore in place multiple sequences such k < i <= j
 * Input:
 *    M - U'dU factorisation (UD format)
 * Output:
 *    M - U'dU recomposition (symmetric)
 */
{
	std::size_t i,j,k;
	const std::size_t n = M.size1();
	assert (n == M.size2());

	// Recompose M = (U'dU) in place
	if (n > 0)
	{
		i = n-1;
		do {
			RowMatrix::Row Mi(M,i);
			// (U' d) row i of lower triangle from upper triangle
			for (j = 0; j < i; ++j)
				Mi[j] = M(j,i) * M(j,j);
			// (U' d) U in place
			j = n-1;
			do { // j>=i
				// Compute matrix product (U'd) row i * U col j
				RowMatrix::value_type Mij = Mi[j];
				if (j > i)					// Optimised handling of 1 in U
					Mij *= Mi[i];
				for (k = 0; k < i; ++k)		// Inner loop k < i <=j, only strict triangular elements
					Mij += Mi[k] * M(k,j);		// M(i,k) element of U'd, M(k,j) element of U
				M(j,i) = Mi[j] = Mij;
			} while (j-- > i);
		} while (i-- > 0);
	}
}


void UdUrecompose (RowMatrix& M)