/* This code is designed to solve min 0.5* s.t. X_ii =1, i=1,2,...,n */ /* and X>=tau*I (symmetric and positive semi-definite) */ /* based on the algorithm in "A Quadratically Convergent Newton Method for */ /* Computing the Nearest Correlation Matrix */ /* By Houduo Qi and Defeng Sun */ /* SIAM J. Matrix Anal. Appl. 28 (2006) 360--385. */ /* This particular C implementation is a result of the summer research project */ /* of Pawel Zaczkowski, Imperial College London, */ /* with Professor Defeng Sun, National University of Singapore */ /* Last modified date: August 11, 2010 */ /* The input argument is the given symmetric G */ /* The outputs are the optimal primal and dual solutions */ /* Diagonal Preconditioner is added */ /* Please send your comments and suggestions to */ /* pawelzaczkowski@cantab.net or matsundf@nus.edu.sg */ /* Warning: Accuracy may not be guaranteed!!!!! */ #include #include "stdio.h" #include "stdlib.h" #include "math.h" #include "time.h" #define PERTURBATION 1.0e-9 /* declare a structure for keeping dimensions and entries of matrices together */ struct matrix { double *entries; int rows; int columns; }; /* function declarations */ void Correlation_Newton(struct matrix* G, struct matrix* X, double* y); /* PURPOSE: calculating the nearest correlation matrix of G */ /* INPUT: struct matrix* G */ /* OUTPUT: struct matrix* X */ /* OUTPUT: double* y */ /* based on the algorithm in "A Quadratically Convergent Newton Method for */ /* Computing the Nearest Correlation Matrix */ /* By Houduo Qi and Defeng Sun */ /* SIAM J. Matrix Anal. Appl. 28 (2006) 360--385. */ void pre_cg(double* b,double tol,int maxit,double* c,struct matrix* Omega12,struct matrix* P, double* p, int* flag, double* relres, int* iterk); /* PURPOSE: PCG method */ /* INPUT: double* b */ /* INPUT: double tol */ /* INPUT: int maxit */ /* INPUT: double* c */ /* INPUT: struct matrix* Omega12 */ /* INPUT: struct matrix* P */ /* OUTPUT: double* p */ /* OUTPUT: int* flag */ /* OUTPUT: double* relres */ /* OUTPUT: int* iterk */ /* This is exactly the algorithm by Hestenes and Stiefel (1952) */ /* An iterative method to solve A(x) =b */ /* The symmetric positive definite matrix M is a preconditioner for A. */ /* See Pages 527 and 534 of Golub and va Loan (1996) */ void precond_matrix (struct matrix* Omega12, struct matrix * P, double* c); /* PURPOSE: generating the diagonal preconditioner */ /* INPUT: struct matrix* Omega12 */ /* INPUT: struct matrix* P */ /* OUTPUT: double* c */ void gradient(double *y, double* lambda, struct matrix* P, double* b0, double* f, double* Fy); /* PURPOSE: generating F(y) */ /* INPUT: double* y */ /* INPUT: double* lambda */ /* INPUT: struct matrix* P */ /* INPUT: double* b0 */ /* OUTPUT: double* f */ /* OUTPUT: double* Fy */ void omega_mat(double* lambda, int n, struct matrix * Omega12); /* PURPOSE: generating the essential part of the first-order difference d */ /* INPUT: double* lambda */ /* INPUT: int n */ /* OUTPUT: struct matrix* Omega12 */ void Jacobian_matrix (double* x, struct matrix* Omega12, struct matrix* P, struct matrix* Ax); /* PURPOSE: generating Jacobian matrix */ /* INPUT: double* x */ /* INPUT: struct matrix* Omega12 */ /* INPUT: struct matrix* P */ /* OUTPUT: struct matrix* Ax */ void printMatrix(struct matrix* M); /* PURPOSE: prints a given matrix on screen */ /* INPUT: struct matrix* M */ double max(double a, double b); double norm(double *a, int n);