-
Lubomir Riha authoredLubomir Riha authored
cs_demo2.c 4.00 KiB
#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) ;
}