dprecond.c

来自「留一模型选择法leave-one-out model selection」· C语言 代码 · 共 35 行

#include <stdlib.h>
#include <stdio.h>

/* LAPACK */
extern int dpotf2_(char *, int *, double *, int *, int *);

double dcholfact(int n, double *A, double *L)
{
	/* if A is p.d. , A = L*L'
	   if A is p.s.d. , A + (1e-3/512/512)*I = L*L'; */  
	int indef, i;
	memcpy(L, A, sizeof(double)*n*n);
	dpotf2_("L", &n, L, &n, &indef);
	if (indef != 0)
	{
		memcpy(L, A, sizeof(double)*n*n);
		for (i=0;i<n;i++)
			L[i*n+i] += (1e-3/512/512); 
		dpotf2_("L", &n, L, &n, &indef);
		if (indef != 0)
			printf("A is not positive semi-definite\n");
		return (1e-3/512/512);
	}
	return 0;
}

double dprecond(int n, double *A, double *C)
{
	/* Given a dense symmetric positive semidefinite matrix A, this
	subroutine computes the precondictioner C. Use full Cholesky
	factorization instead of incomplete Cholesky factorization that
	orginal TRON uses. It is the major difference between the two. */
	return dcholfact(n, A, C);
}