1 /* ========================================================================== */
2 /* === qrtest_C ============================================================= */
3 /* ========================================================================== */
4
5 /* Test the C wrapper functions. */
6
7 #include "SuiteSparseQR_C.h"
8 #define Long SuiteSparse_long
9
10 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
11
qrtest_C(cholmod_sparse * A,double anorm,double errs[5],double maxresid[2][2],cholmod_common * cc)12 void qrtest_C
13 (
14 cholmod_sparse *A,
15 double anorm,
16 double errs [5],
17 double maxresid [2][2],
18 cholmod_common *cc
19 )
20 {
21 cholmod_dense *B, *X, *Resid ;
22 cholmod_sparse *Bsparse, *Xsparse ;
23 SuiteSparseQR_C_factorization *QR ;
24 double resid, one [2] = {1,0}, minusone [2] = {-1,0} ;
25 Long m, n ;
26 #ifndef NEXPERT
27 cholmod_dense *Y ;
28 int split ;
29 #endif
30
31 m = A->nrow ;
32 n = A->ncol ;
33
34 B = cholmod_l_ones (m, 1, A->xtype, cc) ;
35
36 /* X = A\B */
37 X = SuiteSparseQR_C_backslash_default (A, B, cc) ;
38
39 /* Resid = A*X - B */
40 Resid = cholmod_l_copy_dense (B, cc) ;
41 cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ;
42
43 /* resid = norm (Resid,1) */
44 resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ;
45 resid = (resid < 0 || resid != resid) ? 9e99 : resid ;
46 cholmod_l_free_dense (&Resid, cc) ;
47 cholmod_l_free_dense (&X, cc) ;
48
49 maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ;
50 printf ("Resid_C1 %d : %g\n", m>n, resid) ;
51
52 /* X = A\B */
53 Bsparse = cholmod_l_dense_to_sparse (B, 1, cc) ;
54 Xsparse = SuiteSparseQR_C_backslash_sparse (2, -2, A, Bsparse, cc) ;
55 X = cholmod_l_sparse_to_dense (Xsparse, cc) ;
56 cholmod_l_free_sparse (&Bsparse, cc) ;
57 cholmod_l_free_sparse (&Xsparse, cc) ;
58
59 /* Resid = A*X - B */
60 Resid = cholmod_l_copy_dense (B, cc) ;
61 cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ;
62
63 /* resid = norm (Resid,1) */
64 resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ;
65 resid = (resid < 0 || resid != resid) ? 9e99 : resid ;
66 cholmod_l_free_dense (&Resid, cc) ;
67 cholmod_l_free_dense (&X, cc) ;
68
69 maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ;
70 printf ("Resid_C2 %d : %g\n", m>n, resid) ;
71
72 #ifndef NEXPERT
73
74 for (split = 0 ; split <= 1 ; split++)
75 {
76 if (split)
77 {
78 /* split symbolic/numeric QR factorization */
79 QR = SuiteSparseQR_C_symbolic (2, 1, A, cc) ;
80 SuiteSparseQR_C_numeric (-2, A, QR, cc) ;
81 }
82 else
83 {
84 /* QR factorization, single pass */
85 QR = SuiteSparseQR_C_factorize (2, -2, A, cc) ;
86 }
87
88 /* Y = Q'*B */
89 Y = SuiteSparseQR_C_qmult (0, QR, B, cc) ;
90
91 /* X = E*(R\Y) */
92 X = SuiteSparseQR_C_solve (1, QR, Y, cc) ;
93
94 /* Resid = A*X - B */
95 Resid = cholmod_l_copy_dense (B, cc) ;
96 cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ;
97
98 /* resid = norm (Resid,1) */
99 resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ;
100 resid = (resid < 0 || resid != resid) ? 9e99 : resid ;
101 cholmod_l_free_dense (&Resid, cc) ;
102 cholmod_l_free_dense (&X, cc) ;
103
104 maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ;
105 printf ("Resid_C3 %d : %g\n", m>n, resid) ;
106
107 cholmod_l_free_dense (&Y, cc) ;
108 SuiteSparseQR_C_free (&QR, cc) ;
109 }
110 #endif
111
112 cholmod_l_free_dense (&B, cc) ;
113 }
114