Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
amesos_cholmod_complexity.h
Go to the documentation of this file.
1 /* ========================================================================== */
2 /* === Include/cholmod_complexity.h ========================================= */
3 /* ========================================================================== */
4 
5 /* Define operations on pattern, real, complex, and zomplex objects.
6  *
7  * The xtype of an object defines it numerical type. A qttern object has no
8  * numerical values (A->x and A->z are NULL). A real object has no imaginary
9  * qrt (A->x is used, A->z is NULL). A complex object has an imaginary qrt
10  * that is stored interleaved with its real qrt (A->x is of size 2*nz, A->z
11  * is NULL). A zomplex object has both real and imaginary qrts, which are
12  * stored seqrately, as in MATLAB (A->x and A->z are both used).
13  *
14  * XTYPE is CHOLMOD_PATTERN, _REAL, _COMPLEX or _ZOMPLEX, and is the xtype of
15  * the template routine under construction. XTYPE2 is equal to XTYPE, except
16  * if XTYPE is CHOLMOD_PATTERN, in which case XTYPE is CHOLMOD_REAL.
17  * XTYPE and XTYPE2 are defined in cholmod_template.h.
18  */
19 
20 /* -------------------------------------------------------------------------- */
21 /* pattern */
22 /* -------------------------------------------------------------------------- */
23 
24 #define P_TEMPLATE(name) amesos_p_ ## name
25 #define P_ASSIGN2(x,z,p,ax,az,q) x [p] = 1
26 #define P_PRINT(k,x,z,p) PRK(k, ("1"))
27 
28 /* -------------------------------------------------------------------------- */
29 /* real */
30 /* -------------------------------------------------------------------------- */
31 
32 #define R_TEMPLATE(name) amesos_r_ ## name
33 #define R_ASSEMBLE(x,z,p,ax,az,q) x [p] += ax [q]
34 #define R_ASSIGN(x,z,p,ax,az,q) x [p] = ax [q]
35 #define R_ASSIGN_CONJ(x,z,p,ax,az,q) x [p] = ax [q]
36 #define R_ASSIGN_REAL(x,p,ax,q) x [p] = ax [q]
37 #define R_XTYPE_OK(type) ((type) == CHOLMOD_REAL)
38 #define R_IS_NONZERO(ax,az,q) IS_NONZERO (ax [q])
39 #define R_IS_ZERO(ax,az,q) IS_ZERO (ax [q])
40 #define R_IS_ONE(ax,az,q) (ax [q] == 1)
41 #define R_MULT(x,z,p, ax,az,q, bx,bz,r) x [p] = ax [q] * bx [r]
42 #define R_MULTADD(x,z,p, ax,az,q, bx,bz,r) x [p] += ax [q] * bx [r]
43 #define R_MULTSUB(x,z,p, ax,az,q, bx,bz,r) x [p] -= ax [q] * bx [r]
44 #define R_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r) x [p] += ax [q] * bx [r]
45 #define R_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r) x [p] -= ax [q] * bx [r]
46 #define R_ADD(x,z,p, ax,az,q, bx,bz,r) x [p] = ax [q] + bx [r]
47 #define R_ADD_REAL(x,p, ax,q, bx,r) x [p] = ax [q] + bx [r]
48 #define R_CLEAR(x,z,p) x [p] = 0
49 #define R_CLEAR_IMAG(x,z,p)
50 #define R_DIV(x,z,p,ax,az,q) x [p] /= ax [q]
51 #define R_LLDOT(x,p, ax,az,q) x [p] -= ax [q] * ax [q]
52 #define R_PRINT(k,x,z,p) PRK(k, ("%24.16e", x [p]))
53 
54 #define R_DIV_REAL(x,z,p, ax,az,q, bx,r) x [p] = ax [q] / bx [r]
55 #define R_MULT_REAL(x,z,p, ax,az,q, bx,r) x [p] = ax [q] * bx [r]
56 
57 #define R_LDLDOT(x,p, ax,az,q, bx,r) x [p] -=(ax[q] * ax[q])/ bx[r]
58 
59 /* -------------------------------------------------------------------------- */
60 /* complex */
61 /* -------------------------------------------------------------------------- */
62 
63 #define C_TEMPLATE(name) amesos_c_ ## name
64 #define CT_TEMPLATE(name) amesos_ct_ ## name
65 
66 #define C_ASSEMBLE(x,z,p,ax,az,q) \
67  x [2*(p) ] += ax [2*(q) ] ; \
68  x [2*(p)+1] += ax [2*(q)+1]
69 
70 #define C_ASSIGN(x,z,p,ax,az,q) \
71  x [2*(p) ] = ax [2*(q) ] ; \
72  x [2*(p)+1] = ax [2*(q)+1]
73 
74 #define C_ASSIGN_REAL(x,p,ax,q) x [2*(p)] = ax [2*(q)]
75 
76 #define C_ASSIGN_CONJ(x,z,p,ax,az,q) \
77  x [2*(p) ] = ax [2*(q) ] ; \
78  x [2*(p)+1] = -ax [2*(q)+1]
79 
80 #define C_XTYPE_OK(type) ((type) == CHOLMOD_COMPLEX)
81 
82 #define C_IS_NONZERO(ax,az,q) \
83  (IS_NONZERO (ax [2*(q)]) || IS_NONZERO (ax [2*(q)+1]))
84 
85 #define C_IS_ZERO(ax,az,q) \
86  (IS_ZERO (ax [2*(q)]) && IS_ZERO (ax [2*(q)+1]))
87 
88 #define C_IS_ONE(ax,az,q) \
89  ((ax [2*(q)] == 1) && IS_ZERO (ax [2*(q)+1]))
90 
91 #define C_IMAG_IS_NONZERO(ax,az,q) (IS_NONZERO (ax [2*(q)+1]))
92 
93 #define C_MULT(x,z,p, ax,az,q, bx,bz,r) \
94 x [2*(p) ] = ax [2*(q) ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
95 x [2*(p)+1] = ax [2*(q)+1] * bx [2*(r)] + ax [2*(q) ] * bx [2*(r)+1]
96 
97 #define C_MULTADD(x,z,p, ax,az,q, bx,bz,r) \
98 x [2*(p) ] += ax [2*(q) ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
99 x [2*(p)+1] += ax [2*(q)+1] * bx [2*(r)] + ax [2*(q) ] * bx [2*(r)+1]
100 
101 #define C_MULTSUB(x,z,p, ax,az,q, bx,bz,r) \
102 x [2*(p) ] -= ax [2*(q) ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
103 x [2*(p)+1] -= ax [2*(q)+1] * bx [2*(r)] + ax [2*(q) ] * bx [2*(r)+1]
104 
105 /* s += conj(a)*b */
106 #define C_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r) \
107 x [2*(p) ] += ax [2*(q) ] * bx [2*(r)] + ax [2*(q)+1] * bx [2*(r)+1] ; \
108 x [2*(p)+1] += (-ax [2*(q)+1]) * bx [2*(r)] + ax [2*(q) ] * bx [2*(r)+1]
109 
110 /* s -= conj(a)*b */
111 #define C_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r) \
112 x [2*(p) ] -= ax [2*(q) ] * bx [2*(r)] + ax [2*(q)+1] * bx [2*(r)+1] ; \
113 x [2*(p)+1] -= (-ax [2*(q)+1]) * bx [2*(r)] + ax [2*(q) ] * bx [2*(r)+1]
114 
115 #define C_ADD(x,z,p, ax,az,q, bx,bz,r) \
116  x [2*(p) ] = ax [2*(q) ] + bx [2*(r) ] ; \
117  x [2*(p)+1] = ax [2*(q)+1] + bx [2*(r)+1]
118 
119 #define C_ADD_REAL(x,p, ax,q, bx,r) \
120  x [2*(p)] = ax [2*(q)] + bx [2*(r)]
121 
122 #define C_CLEAR(x,z,p) \
123  x [2*(p) ] = 0 ; \
124  x [2*(p)+1] = 0
125 
126 #define C_CLEAR_IMAG(x,z,p) \
127  x [2*(p)+1] = 0
128 
129 /* s = s / a */
130 #define C_DIV(x,z,p,ax,az,q) \
131  Common->complex_divide ( \
132  x [2*(p)], x [2*(p)+1], \
133  ax [2*(q)], ax [2*(q)+1], \
134  &x [2*(p)], &x [2*(p)+1])
135 
136 /* s -= conj(a)*a ; note that the result of conj(a)*a is real */
137 #define C_LLDOT(x,p, ax,az,q) \
138  x [2*(p)] -= ax [2*(q)] * ax [2*(q)] + ax [2*(q)+1] * ax [2*(q)+1]
139 
140 #define C_PRINT(k,x,z,p) PRK(k, ("(%24.16e,%24.16e)", x [2*(p)], x [2*(p)+1]))
141 
142 #define C_DIV_REAL(x,z,p, ax,az,q, bx,r) \
143  x [2*(p) ] = ax [2*(q) ] / bx [2*(r)] ; \
144  x [2*(p)+1] = ax [2*(q)+1] / bx [2*(r)]
145 
146 #define C_MULT_REAL(x,z,p, ax,az,q, bx,r) \
147  x [2*(p) ] = ax [2*(q) ] * bx [2*(r)] ; \
148  x [2*(p)+1] = ax [2*(q)+1] * bx [2*(r)]
149 
150 /* s -= conj(a)*a/t */
151 #define C_LDLDOT(x,p, ax,az,q, bx,r) \
152  x [2*(p)] -= (ax [2*(q)] * ax [2*(q)] + ax [2*(q)+1] * ax [2*(q)+1]) / bx[r]
153 
154 /* -------------------------------------------------------------------------- */
155 /* zomplex */
156 /* -------------------------------------------------------------------------- */
157 
158 #define Z_TEMPLATE(name) amesos_z_ ## name
159 #define ZT_TEMPLATE(name) amesos_zt_ ## name
160 
161 #define Z_ASSEMBLE(x,z,p,ax,az,q) \
162  x [p] += ax [q] ; \
163  z [p] += az [q]
164 
165 #define Z_ASSIGN(x,z,p,ax,az,q) \
166  x [p] = ax [q] ; \
167  z [p] = az [q]
168 
169 #define Z_ASSIGN_REAL(x,p,ax,q) x [p] = ax [q]
170 
171 #define Z_ASSIGN_CONJ(x,z,p,ax,az,q) \
172  x [p] = ax [q] ; \
173  z [p] = -az [q]
174 
175 #define Z_XTYPE_OK(type) ((type) == CHOLMOD_ZOMPLEX)
176 
177 #define Z_IS_NONZERO(ax,az,q) \
178  (IS_NONZERO (ax [q]) || IS_NONZERO (az [q]))
179 
180 #define Z_IS_ZERO(ax,az,q) \
181  (IS_ZERO (ax [q]) && IS_ZERO (az [q]))
182 
183 #define Z_IS_ONE(ax,az,q) \
184  ((ax [q] == 1) && IS_ZERO (az [q]))
185 
186 #define Z_IMAG_IS_NONZERO(ax,az,q) (IS_NONZERO (az [q]))
187 
188 #define Z_MULT(x,z,p, ax,az,q, bx,bz,r) \
189  x [p] = ax [q] * bx [r] - az [q] * bz [r] ; \
190  z [p] = az [q] * bx [r] + ax [q] * bz [r]
191 
192 #define Z_MULTADD(x,z,p, ax,az,q, bx,bz,r) \
193  x [p] += ax [q] * bx [r] - az [q] * bz [r] ; \
194  z [p] += az [q] * bx [r] + ax [q] * bz [r]
195 
196 #define Z_MULTSUB(x,z,p, ax,az,q, bx,bz,r) \
197  x [p] -= ax [q] * bx [r] - az [q] * bz [r] ; \
198  z [p] -= az [q] * bx [r] + ax [q] * bz [r]
199 
200 #define Z_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r) \
201  x [p] += ax [q] * bx [r] + az [q] * bz [r] ; \
202  z [p] += (-az [q]) * bx [r] + ax [q] * bz [r]
203 
204 #define Z_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r) \
205  x [p] -= ax [q] * bx [r] + az [q] * bz [r] ; \
206  z [p] -= (-az [q]) * bx [r] + ax [q] * bz [r]
207 
208 #define Z_ADD(x,z,p, ax,az,q, bx,bz,r) \
209  x [p] = ax [q] + bx [r] ; \
210  z [p] = az [q] + bz [r]
211 
212 #define Z_ADD_REAL(x,p, ax,q, bx,r) \
213  x [p] = ax [q] + bx [r]
214 
215 #define Z_CLEAR(x,z,p) \
216  x [p] = 0 ; \
217  z [p] = 0
218 
219 #define Z_CLEAR_IMAG(x,z,p) \
220  z [p] = 0
221 
222 /* s = s/a */
223 #define Z_DIV(x,z,p,ax,az,q) \
224  Common->complex_divide (x [p], z [p], ax [q], az [q], &x [p], &z [p])
225 
226 /* s -= conj(a)*a ; note that the result of conj(a)*a is real */
227 #define Z_LLDOT(x,p, ax,az,q) \
228  x [p] -= ax [q] * ax [q] + az [q] * az [q]
229 
230 #define Z_PRINT(k,x,z,p) PRK(k, ("(%24.16e,%24.16e)", x [p], z [p]))
231 
232 #define Z_DIV_REAL(x,z,p, ax,az,q, bx,r) \
233  x [p] = ax [q] / bx [r] ; \
234  z [p] = az [q] / bx [r]
235 
236 #define Z_MULT_REAL(x,z,p, ax,az,q, bx,r) \
237  x [p] = ax [q] * bx [r] ; \
238  z [p] = az [q] * bx [r]
239 
240 /* s -= conj(a)*a/t */
241 #define Z_LDLDOT(x,p, ax,az,q, bx,r) \
242  x [p] -= (ax [q] * ax [q] + az [q] * az [q]) / bx[r]
243 
244 /* -------------------------------------------------------------------------- */
245 /* all classes */
246 /* -------------------------------------------------------------------------- */
247 
248 /* Check if A->xtype and the two arrays A->x and A->z are valid. Set status to
249  * invalid, unless status is already "out of memory". A can be a sparse matrix,
250  * dense matrix, factor, or triplet. */
251 
252 #define RETURN_IF_XTYPE_INVALID(A,xtype1,xtype2,result) \
253 { \
254  if ((A)->xtype < (xtype1) || (A)->xtype > (xtype2) || \
255  ((A)->xtype != CHOLMOD_PATTERN && ((A)->x) == NULL) || \
256  ((A)->xtype == CHOLMOD_ZOMPLEX && ((A)->z) == NULL)) \
257  { \
258  if (Common->status != CHOLMOD_OUT_OF_MEMORY) \
259  { \
260  ERROR (CHOLMOD_INVALID, "invalid xtype") ; \
261  } \
262  return (result) ; \
263  } \
264 }