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) 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) 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) c_ ## name 64 #define CT_TEMPLATE(name) 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 SuiteSparse_config.divcomplex_func ( \ 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) z_ ## name 159 #define ZT_TEMPLATE(name) 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 SuiteSparse_config.divcomplex_func \ 225 (x [p], z [p], ax [q], az [q], &x [p], &z [p]) 226 227 /* s -= conj(a)*a ; note that the result of conj(a)*a is real */ 228 #define Z_LLDOT(x,p, ax,az,q) \ 229 x [p] -= ax [q] * ax [q] + az [q] * az [q] 230 231 #define Z_PRINT(k,x,z,p) PRK(k, ("(%24.16e,%24.16e)", x [p], z [p])) 232 233 #define Z_DIV_REAL(x,z,p, ax,az,q, bx,r) \ 234 x [p] = ax [q] / bx [r] ; \ 235 z [p] = az [q] / bx [r] 236 237 #define Z_MULT_REAL(x,z,p, ax,az,q, bx,r) \ 238 x [p] = ax [q] * bx [r] ; \ 239 z [p] = az [q] * bx [r] 240 241 /* s -= conj(a)*a/t */ 242 #define Z_LDLDOT(x,p, ax,az,q, bx,r) \ 243 x [p] -= (ax [q] * ax [q] + az [q] * az [q]) / bx[r] 244 245 /* -------------------------------------------------------------------------- */ 246 /* all classes */ 247 /* -------------------------------------------------------------------------- */ 248 249 /* Check if A->xtype and the two arrays A->x and A->z are valid. Set status to 250 * invalid, unless status is already "out of memory". A can be a sparse matrix, 251 * dense matrix, factor, or triplet. */ 252 253 #define RETURN_IF_XTYPE_INVALID(A,xtype1,xtype2,result) \ 254 { \ 255 if ((A)->xtype < (xtype1) || (A)->xtype > (xtype2) || \ 256 ((A)->xtype != CHOLMOD_PATTERN && ((A)->x) == NULL) || \ 257 ((A)->xtype == CHOLMOD_ZOMPLEX && ((A)->z) == NULL)) \ 258 { \ 259 if (Common->status != CHOLMOD_OUT_OF_MEMORY) \ 260 { \ 261 ERROR (CHOLMOD_INVALID, "invalid xtype") ; \ 262 } \ 263 return (result) ; \ 264 } \ 265 } 266