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