suitesparse: cholmod_complexity.h Source File

00001 /* ========================================================================== */
00002 /* === Include/cholmod_complexity.h ========================================= */
00003 /* ========================================================================== */
00004 
00005 /* Define operations on pattern, real, complex, and zomplex objects.
00006  *
00007  * The xtype of an object defines it numerical type.  A qttern object has no
00008  * numerical values (A->x and A->z are NULL).  A real object has no imaginary
00009  * qrt (A->x is used, A->z is NULL).  A complex object has an imaginary qrt
00010  * that is stored interleaved with its real qrt (A->x is of size 2*nz, A->z
00011  * is NULL).  A zomplex object has both real and imaginary qrts, which are
00012  * stored seqrately, as in MATLAB (A->x and A->z are both used).
00013  *
00014  * XTYPE is CHOLMOD_PATTERN, _REAL, _COMPLEX or _ZOMPLEX, and is the xtype of
00015  * the template routine under construction.  XTYPE2 is equal to XTYPE, except
00016  * if XTYPE is CHOLMOD_PATTERN, in which case XTYPE is CHOLMOD_REAL.
00017  * XTYPE and XTYPE2 are defined in cholmod_template.h.  
00018  */
00019 
00020 /* -------------------------------------------------------------------------- */
00021 /* pattern */
00022 /* -------------------------------------------------------------------------- */
00023 
00024 #define P_TEMPLATE(name)                p_ ## name
00025 #define P_ASSIGN2(x,z,p,ax,az,q)        x [p] = 1
00026 #define P_PRINT(k,x,z,p)                PRK(k, ("1"))
00027 
00028 /* -------------------------------------------------------------------------- */
00029 /* real */
00030 /* -------------------------------------------------------------------------- */
00031 
00032 #define R_TEMPLATE(name)                        r_ ## name
00033 #define R_ASSEMBLE(x,z,p,ax,az,q)               x [p] += ax [q]
00034 #define R_ASSIGN(x,z,p,ax,az,q)                 x [p]  = ax [q]
00035 #define R_ASSIGN_CONJ(x,z,p,ax,az,q)            x [p]  = ax [q]
00036 #define R_ASSIGN_REAL(x,p,ax,q)                 x [p]  = ax [q]
00037 #define R_XTYPE_OK(type)                        ((type) == CHOLMOD_REAL)
00038 #define R_IS_NONZERO(ax,az,q)                   IS_NONZERO (ax [q])
00039 #define R_IS_ZERO(ax,az,q)                      IS_ZERO (ax [q])
00040 #define R_IS_ONE(ax,az,q)                       (ax [q] == 1)
00041 #define R_MULT(x,z,p, ax,az,q, bx,bz,r)         x [p]  = ax [q] * bx [r]
00042 #define R_MULTADD(x,z,p, ax,az,q, bx,bz,r)      x [p] += ax [q] * bx [r]
00043 #define R_MULTSUB(x,z,p, ax,az,q, bx,bz,r)      x [p] -= ax [q] * bx [r]
00044 #define R_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r)  x [p] += ax [q] * bx [r]
00045 #define R_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r)  x [p] -= ax [q] * bx [r]
00046 #define R_ADD(x,z,p, ax,az,q, bx,bz,r)          x [p]  = ax [q] + bx [r]
00047 #define R_ADD_REAL(x,p, ax,q, bx,r)             x [p]  = ax [q] + bx [r]
00048 #define R_CLEAR(x,z,p)                          x [p]  = 0
00049 #define R_CLEAR_IMAG(x,z,p)
00050 #define R_DIV(x,z,p,ax,az,q)                    x [p] /= ax [q]
00051 #define R_LLDOT(x,p, ax,az,q)                   x [p] -= ax [q] * ax [q]
00052 #define R_PRINT(k,x,z,p)                        PRK(k, ("%24.16e", x [p]))
00053 
00054 #define R_DIV_REAL(x,z,p, ax,az,q, bx,r)        x [p] = ax [q] / bx [r]
00055 #define R_MULT_REAL(x,z,p, ax,az,q, bx,r)       x [p] = ax [q] * bx [r]
00056 
00057 #define R_LDLDOT(x,p, ax,az,q, bx,r)            x [p] -=(ax[q] * ax[q])/ bx[r]
00058 
00059 /* -------------------------------------------------------------------------- */
00060 /* complex */
00061 /* -------------------------------------------------------------------------- */
00062 
00063 #define C_TEMPLATE(name)                c_ ## name
00064 #define CT_TEMPLATE(name)               ct_ ## name
00065 
00066 #define C_ASSEMBLE(x,z,p,ax,az,q) \
00067     x [2*(p)  ] += ax [2*(q)  ] ; \
00068     x [2*(p)+1] += ax [2*(q)+1]
00069 
00070 #define C_ASSIGN(x,z,p,ax,az,q) \
00071     x [2*(p)  ] = ax [2*(q)  ] ; \
00072     x [2*(p)+1] = ax [2*(q)+1]
00073 
00074 #define C_ASSIGN_REAL(x,p,ax,q)                 x [2*(p)]  = ax [2*(q)]
00075 
00076 #define C_ASSIGN_CONJ(x,z,p,ax,az,q) \
00077     x [2*(p)  ] =  ax [2*(q)  ] ; \
00078     x [2*(p)+1] = -ax [2*(q)+1]
00079 
00080 #define C_XTYPE_OK(type)                ((type) == CHOLMOD_COMPLEX)
00081 
00082 #define C_IS_NONZERO(ax,az,q) \
00083     (IS_NONZERO (ax [2*(q)]) || IS_NONZERO (ax [2*(q)+1]))
00084 
00085 #define C_IS_ZERO(ax,az,q) \
00086     (IS_ZERO (ax [2*(q)]) && IS_ZERO (ax [2*(q)+1]))
00087 
00088 #define C_IS_ONE(ax,az,q) \
00089     ((ax [2*(q)] == 1) && IS_ZERO (ax [2*(q)+1]))
00090 
00091 #define C_IMAG_IS_NONZERO(ax,az,q)  (IS_NONZERO (ax [2*(q)+1]))
00092 
00093 #define C_MULT(x,z,p, ax,az,q, bx,bz,r) \
00094 x [2*(p)  ] = ax [2*(q)  ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
00095 x [2*(p)+1] = ax [2*(q)+1] * bx [2*(r)] + ax [2*(q)  ] * bx [2*(r)+1]
00096 
00097 #define C_MULTADD(x,z,p, ax,az,q, bx,bz,r) \
00098 x [2*(p)  ] += ax [2*(q)  ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
00099 x [2*(p)+1] += ax [2*(q)+1] * bx [2*(r)] + ax [2*(q)  ] * bx [2*(r)+1]
00100 
00101 #define C_MULTSUB(x,z,p, ax,az,q, bx,bz,r) \
00102 x [2*(p)  ] -= ax [2*(q)  ] * bx [2*(r)] - ax [2*(q)+1] * bx [2*(r)+1] ; \
00103 x [2*(p)+1] -= ax [2*(q)+1] * bx [2*(r)] + ax [2*(q)  ] * bx [2*(r)+1]
00104 
00105 /* s += conj(a)*b */
00106 #define C_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r) \
00107 x [2*(p)  ] +=   ax [2*(q)  ]  * bx [2*(r)] + ax [2*(q)+1] * bx [2*(r)+1] ; \
00108 x [2*(p)+1] += (-ax [2*(q)+1]) * bx [2*(r)] + ax [2*(q)  ] * bx [2*(r)+1]
00109 
00110 /* s -= conj(a)*b */
00111 #define C_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r) \
00112 x [2*(p)  ] -=   ax [2*(q)  ]  * bx [2*(r)] + ax [2*(q)+1] * bx [2*(r)+1] ; \
00113 x [2*(p)+1] -= (-ax [2*(q)+1]) * bx [2*(r)] + ax [2*(q)  ] * bx [2*(r)+1]
00114 
00115 #define C_ADD(x,z,p, ax,az,q, bx,bz,r) \
00116     x [2*(p)  ] = ax [2*(q)  ] + bx [2*(r)  ] ; \
00117     x [2*(p)+1] = ax [2*(q)+1] + bx [2*(r)+1]
00118 
00119 #define C_ADD_REAL(x,p, ax,q, bx,r) \
00120     x [2*(p)] = ax [2*(q)] + bx [2*(r)]
00121 
00122 #define C_CLEAR(x,z,p) \
00123     x [2*(p)  ] = 0 ; \
00124     x [2*(p)+1] = 0
00125 
00126 #define C_CLEAR_IMAG(x,z,p) \
00127     x [2*(p)+1] = 0
00128 
00129 /* s = s / a */
00130 #define C_DIV(x,z,p,ax,az,q) \
00131     Common->complex_divide ( \
00132               x [2*(p)],  x [2*(p)+1], \
00133              ax [2*(q)], ax [2*(q)+1], \
00134              &x [2*(p)], &x [2*(p)+1])
00135 
00136 /* s -= conj(a)*a ; note that the result of conj(a)*a is real */
00137 #define C_LLDOT(x,p, ax,az,q) \
00138     x [2*(p)] -= ax [2*(q)] * ax [2*(q)] + ax [2*(q)+1] * ax [2*(q)+1]
00139 
00140 #define C_PRINT(k,x,z,p) PRK(k, ("(%24.16e,%24.16e)", x [2*(p)], x [2*(p)+1]))
00141 
00142 #define C_DIV_REAL(x,z,p, ax,az,q, bx,r) \
00143     x [2*(p)  ] = ax [2*(q)  ] / bx [2*(r)] ; \
00144     x [2*(p)+1] = ax [2*(q)+1] / bx [2*(r)]
00145 
00146 #define C_MULT_REAL(x,z,p, ax,az,q, bx,r) \
00147     x [2*(p)  ] = ax [2*(q)  ] * bx [2*(r)] ; \
00148     x [2*(p)+1] = ax [2*(q)+1] * bx [2*(r)]
00149 
00150 /* s -= conj(a)*a/t */
00151 #define C_LDLDOT(x,p, ax,az,q, bx,r) \
00152     x [2*(p)] -= (ax [2*(q)] * ax [2*(q)] + ax [2*(q)+1] * ax [2*(q)+1]) / bx[r]
00153 
00154 /* -------------------------------------------------------------------------- */
00155 /* zomplex */
00156 /* -------------------------------------------------------------------------- */
00157 
00158 #define Z_TEMPLATE(name)                z_ ## name
00159 #define ZT_TEMPLATE(name)               zt_ ## name
00160 
00161 #define Z_ASSEMBLE(x,z,p,ax,az,q) \
00162     x [p] += ax [q] ; \
00163     z [p] += az [q]
00164 
00165 #define Z_ASSIGN(x,z,p,ax,az,q) \
00166     x [p] = ax [q] ; \
00167     z [p] = az [q]
00168 
00169 #define Z_ASSIGN_REAL(x,p,ax,q)                 x [p]  = ax [q]
00170 
00171 #define Z_ASSIGN_CONJ(x,z,p,ax,az,q) \
00172     x [p] =  ax [q] ; \
00173     z [p] = -az [q]
00174 
00175 #define Z_XTYPE_OK(type)                ((type) == CHOLMOD_ZOMPLEX)
00176 
00177 #define Z_IS_NONZERO(ax,az,q) \
00178     (IS_NONZERO (ax [q]) || IS_NONZERO (az [q]))
00179 
00180 #define Z_IS_ZERO(ax,az,q) \
00181     (IS_ZERO (ax [q]) && IS_ZERO (az [q]))
00182 
00183 #define Z_IS_ONE(ax,az,q) \
00184     ((ax [q] == 1) && IS_ZERO (az [q]))
00185 
00186 #define Z_IMAG_IS_NONZERO(ax,az,q)  (IS_NONZERO (az [q]))
00187 
00188 #define Z_MULT(x,z,p, ax,az,q, bx,bz,r) \
00189     x [p] = ax [q] * bx [r] - az [q] * bz [r] ; \
00190     z [p] = az [q] * bx [r] + ax [q] * bz [r]
00191 
00192 #define Z_MULTADD(x,z,p, ax,az,q, bx,bz,r) \
00193     x [p] += ax [q] * bx [r] - az [q] * bz [r] ; \
00194     z [p] += az [q] * bx [r] + ax [q] * bz [r]
00195 
00196 #define Z_MULTSUB(x,z,p, ax,az,q, bx,bz,r) \
00197     x [p] -= ax [q] * bx [r] - az [q] * bz [r] ; \
00198     z [p] -= az [q] * bx [r] + ax [q] * bz [r]
00199 
00200 #define Z_MULTADDCONJ(x,z,p, ax,az,q, bx,bz,r) \
00201     x [p] +=   ax [q]  * bx [r] + az [q] * bz [r] ; \
00202     z [p] += (-az [q]) * bx [r] + ax [q] * bz [r]
00203 
00204 #define Z_MULTSUBCONJ(x,z,p, ax,az,q, bx,bz,r) \
00205     x [p] -=   ax [q]  * bx [r] + az [q] * bz [r] ; \
00206     z [p] -= (-az [q]) * bx [r] + ax [q] * bz [r]
00207 
00208 #define Z_ADD(x,z,p, ax,az,q, bx,bz,r) \
00209         x [p] = ax [q] + bx [r] ; \
00210         z [p] = az [q] + bz [r]
00211 
00212 #define Z_ADD_REAL(x,p, ax,q, bx,r) \
00213         x [p] = ax [q] + bx [r]
00214 
00215 #define Z_CLEAR(x,z,p) \
00216     x [p] = 0 ; \
00217     z [p] = 0
00218 
00219 #define Z_CLEAR_IMAG(x,z,p) \
00220     z [p] = 0
00221 
00222 /* s = s/a */
00223 #define Z_DIV(x,z,p,ax,az,q) \
00224     Common->complex_divide (x [p], z [p], ax [q], az [q], &x [p], &z [p])
00225 
00226 /* s -= conj(a)*a ; note that the result of conj(a)*a is real */
00227 #define Z_LLDOT(x,p, ax,az,q) \
00228     x [p] -= ax [q] * ax [q] + az [q] * az [q]
00229 
00230 #define Z_PRINT(k,x,z,p)        PRK(k, ("(%24.16e,%24.16e)", x [p], z [p]))
00231 
00232 #define Z_DIV_REAL(x,z,p, ax,az,q, bx,r) \
00233     x [p] = ax [q] / bx [r] ; \
00234     z [p] = az [q] / bx [r]
00235 
00236 #define Z_MULT_REAL(x,z,p, ax,az,q, bx,r) \
00237     x [p] = ax [q] * bx [r] ; \
00238     z [p] = az [q] * bx [r]
00239 
00240 /* s -= conj(a)*a/t */
00241 #define Z_LDLDOT(x,p, ax,az,q, bx,r) \
00242     x [p] -= (ax [q] * ax [q] + az [q] * az [q]) / bx[r]
00243 
00244 /* -------------------------------------------------------------------------- */
00245 /* all classes */
00246 /* -------------------------------------------------------------------------- */
00247 
00248 /* Check if A->xtype and the two arrays A->x and A->z are valid.  Set status to
00249  * invalid, unless status is already "out of memory".  A can be a sparse matrix,
00250  * dense matrix, factor, or triplet. */
00251 
00252 #define RETURN_IF_XTYPE_INVALID(A,xtype1,xtype2,result) \
00253 { \
00254     if ((A)->xtype < (xtype1) || (A)->xtype > (xtype2) || \
00255         ((A)->xtype != CHOLMOD_PATTERN && ((A)->x) == NULL) || \
00256         ((A)->xtype == CHOLMOD_ZOMPLEX && ((A)->z) == NULL)) \
00257     { \
00258         if (Common->status != CHOLMOD_OUT_OF_MEMORY) \
00259         { \
00260             ERROR (CHOLMOD_INVALID, "invalid xtype") ; \
00261         } \
00262         return (result) ; \
00263     } \
00264 }