c_cpp 计算Cholesky分解

Posted 2021-05-14

void calcCholesky( int nP, double * fCov, double* fCovSqrt ){

		double *C = fCovSqrt;
		double *V = fCov;

		// calculate sqrt(V) as lower diagonal matrix
		for( int i = 0; i < nP; ++i ) {
			for( int j = 0; j < nP; ++j ) {
				C[i*nP+j] = 0;
			}
		}

		for( int j = 0; j < nP; ++j ) {
			// diagonal terms first
			double Ck = 0;
			for( int k = 0; k < j; ++k ) {
				Ck += C[j*nP+k] * C[j*nP+k];
			} // k
			C[j*nP+j] = sqrt( fabs( V[j*nP+j] - Ck ) );

			// off-diagonal terms
			for( int i = j+1; i < nP; ++i ) {
				Ck = 0;
				for( int k = 0; k < j; ++k ) {
					Ck += C[i*nP+k] * C[j*nP+k];
				} //k
				if( C[ j * nP + j ] != 0 ) 
					C[ i * nP + j ] = ( V[ i * nP + j ] - Ck ) / C[ j * nP + j ];
				else 
					C[ i * nP + j ] = 0;
			}// i
		} // j
	} // calcCholesky