#include "cs_demo.h" #include <time.h> /* cs_demo2: read a matrix and solve a linear system */ static double tic (void) { return (clock () / (double) CLOCKS_PER_SEC) ; } static double toc (double t) { double s = tic () ; return (CS_MAX (0, s-t)) ; } static void print_order (csi order) { switch (order) { case 0: printf ("natural ") ; break ; case 1: printf ("amd(A+A') ") ; break ; case 2: printf ("amd(S'*S) ") ; break ; case 3: printf ("amd(A'*A) ") ; break ; } } static void rhs (double *x, double *b, csi m) { csi i ; for (i = 0 ; i < m ; i++) b [i] = 1 + ((double) i) / m ; for (i = 0 ; i < m ; i++) x [i] = b [i] ; } static double norm (double *x, csi n) { csi i ; double normx = 0 ; for (i = 0 ; i < n ; i++) normx = CS_MAX (normx, fabs (x [i])) ; return (normx) ; } static void print_resid (csi ok, cs *A, double *x, double *b, double *resid) { csi i, m, n ; if (!ok) { printf (" (failed)\n") ; return ; } m = A->m ; n = A->n ; for (i = 0 ; i < m ; i++) resid [i] = -b [i] ; /* resid = -b */ cs_gaxpy (A, x, resid) ; /* resid = resid + A*x */ printf ("resid: %8.2e\n", norm (resid,m) / ((n == 0) ? 1 : (cs_norm (A) * norm (x,n) + norm (b,m)))) ; } int main (void) { // problem *Prob = get_problem (stdin, 1e-14) ; // demo2 (Prob) ; // free_problem (Prob) ; problem *Prob = get_problem (stdin, 1e-14) ; cs *A, *C ; double *b, *x, *resid, t, tol ; csi k, m, n, ok, order, nb, ns, *r, *s, *rr, sprank ; csd *D ; if (!Prob) return (0) ; A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid; m = A->m ; n = A->n ; tol = Prob->sym ? 0.001 : 1 ; /* partial pivoting tolerance */ D = cs_dmperm (C, 1) ; /* randomized dmperm analysis */ if (!D) return (0) ; nb = D->nb ; r = D->r ; s = D->s ; rr = D->rr ; sprank = rr [3] ; for (ns = 0, k = 0 ; k < nb ; k++) { ns += ((r [k+1] == r [k]+1) && (s [k+1] == s [k]+1)) ; } printf ("blocks: %g singletons: %g structural rank: %g\n", (double) nb, (double) ns, (double) sprank) ; cs_dfree (D) ; // for (order = 0 ; order <= 3 ; order += 3) /* natural and amd(A'*A) */ // { // if (!order && m > 1000) continue ; // printf ("QR ") ; // print_order (order) ; // rhs (x, b, m) ; /* compute right-hand side */ // t = tic () ; // ok = cs_qrsol (order, C, x) ; // min norm(Ax-b) with QR // printf ("time: %8.2f ", toc (t)) ; // print_resid (ok, C, x, b, resid) ; /* print residual */ // } // if (m != n || sprank < n) return (1) ; // return if rect. or singular // for (order = 0 ; order <= 3 ; order++) /* try all orderings */ // { // if (!order && m > 1000) continue ; // printf ("LU ") ; // print_order (order) ; // rhs (x, b, m) ; /* compute right-hand side */ // t = tic () ; // ok = cs_lusol (order, C, x, tol) ; /* solve Ax=b with LU */ // printf ("time: %8.2f ", toc (t)) ; // print_resid (ok, C, x, b, resid) ; /* print residual */ // } if (!Prob->sym) return (1) ; for (order = 0 ; order <= 1 ; order++) /* natural and amd(A+A') */ { if (!order && m > 1000) continue ; printf ("Chol \n") ; print_order (order) ; rhs (x, b, m) ; /* compute right-hand side */ t = tic () ; ok = cs_cholsol (order, C, x) ; /* solve Ax=b with Cholesky */ printf ("time: %8.2f ", toc (t)) ; print_resid (ok, C, x, b, resid) ; /* print residual */ } // return (1) ; // cs *T, *A, *Eye, *AT, *C, *D ; // csi i, m ; // T = cs_load2 (stdin) ; load triplet matrix T from stdin // printf ("T:\n") ; cs_print (T, 0) ; /* print T */ return (0) ; }