Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_version.h
Go to the documentation of this file.
00001 #ifndef AMESOS_KLU_VERSION_H
00002 #define AMESOS_KLU_VERSION_H
00003 
00004 #ifdef DLONG
00005 #define Int UF_long
00006 #define Int_id UF_long_id
00007 #define Int_MAX UF_long_max
00008 #else
00009 #define Int int
00010 #define Int_id "%d"
00011 #define Int_MAX INT_MAX
00012 #endif
00013 
00014 #define NPRINT  
00015 
00016 #define BYTES(type,n) (sizeof (type) * (n))
00017 #define CEILING(b,u)  (((b)+(u)-1) / (u))
00018 #define UNITS(type,n) (CEILING (BYTES (type,n), sizeof (Unit)))
00019 #define DUNITS(type,n) (ceil (BYTES (type, (double) n) / sizeof (Unit)))
00020 
00021 #define GET_I_POINTER(LU, Xip, Xi, k) \
00022 { \
00023     Xi = (Int *) (LU + Xip [k]) ; \
00024 }
00025 
00026 #define GET_X_POINTER(LU, Xip, Xlen, Xx, k) \
00027 { \
00028     Xx = (Entry *) (LU + Xip [k] + UNITS (Int, Xlen [k])) ; \
00029 }
00030 
00031 #define GET_POINTER(LU, Xip, Xlen, Xi, Xx, k, xlen) \
00032 { \
00033     Unit *xp = LU + Xip [k] ; \
00034     xlen = Xlen [k] ; \
00035     Xi = (Int *) xp ; \
00036     Xx = (Entry *) (xp + UNITS (Int, xlen)) ; \
00037 }
00038 
00039 /* function names */
00040 #ifdef COMPLEX 
00041 
00042 #ifdef DLONG
00043 
00044 #define KLU_scale amesos_klu_zl_scale
00045 #define KLU_solve amesos_klu_zl_solve
00046 #define KLU_tsolve amesos_klu_zl_tsolve
00047 #define KLU_free_numeric amesos_klu_zl_free_numeric
00048 #define KLU_factor amesos_klu_zl_factor
00049 #define KLU_refactor amesos_klu_zl_refactor
00050 #define KLU_kernel_factor amesos_klu_zl_kernel_factor 
00051 #define KLU_lsolve amesos_klu_zl_lsolve
00052 #define KLU_ltsolve amesos_klu_zl_ltsolve
00053 #define KLU_usolve amesos_klu_zl_usolve
00054 #define KLU_utsolve amesos_klu_zl_utsolve
00055 #define KLU_kernel amesos_klu_zl_kernel
00056 #define KLU_valid amesos_klu_zl_valid
00057 #define KLU_valid_LU amesos_klu_zl_valid_LU
00058 #define KLU_sort amesos_klu_zl_sort
00059 #define KLU_rgrowth amesos_klu_zl_rgrowth
00060 #define KLU_rcond amesos_klu_zl_rcond
00061 #define KLU_extract amesos_klu_zl_extract
00062 #define KLU_condest amesos_klu_zl_condest
00063 #define KLU_flops amesos_klu_zl_flops
00064 
00065 #else
00066 
00067 #define KLU_scale amesos_klu_z_scale
00068 #define KLU_solve amesos_klu_z_solve
00069 #define KLU_tsolve amesos_klu_z_tsolve
00070 #define KLU_free_numeric amesos_klu_z_free_numeric
00071 #define KLU_factor amesos_klu_z_factor
00072 #define KLU_refactor amesos_klu_z_refactor
00073 #define KLU_kernel_factor amesos_klu_z_kernel_factor 
00074 #define KLU_lsolve amesos_klu_z_lsolve
00075 #define KLU_ltsolve amesos_klu_z_ltsolve
00076 #define KLU_usolve amesos_klu_z_usolve
00077 #define KLU_utsolve amesos_klu_z_utsolve
00078 #define KLU_kernel amesos_klu_z_kernel
00079 #define KLU_valid amesos_klu_z_valid
00080 #define KLU_valid_LU amesos_klu_z_valid_LU
00081 #define KLU_sort amesos_klu_z_sort
00082 #define KLU_rgrowth amesos_klu_z_rgrowth
00083 #define KLU_rcond amesos_klu_z_rcond
00084 #define KLU_extract amesos_klu_z_extract
00085 #define KLU_condest amesos_klu_z_condest
00086 #define KLU_flops amesos_klu_z_flops
00087 
00088 #endif
00089 
00090 #else
00091 
00092 #ifdef DLONG
00093 
00094 #define KLU_scale amesos_klu_l_scale
00095 #define KLU_solve amesos_klu_l_solve
00096 #define KLU_tsolve amesos_klu_l_tsolve
00097 #define KLU_free_numeric amesos_klu_l_free_numeric
00098 #define KLU_factor amesos_klu_l_factor
00099 #define KLU_refactor amesos_klu_l_refactor
00100 #define KLU_kernel_factor amesos_klu_l_kernel_factor 
00101 #define KLU_lsolve amesos_klu_l_lsolve
00102 #define KLU_ltsolve amesos_klu_l_ltsolve
00103 #define KLU_usolve amesos_klu_l_usolve
00104 #define KLU_utsolve amesos_klu_l_utsolve
00105 #define KLU_kernel amesos_klu_l_kernel
00106 #define KLU_valid amesos_klu_l_valid
00107 #define KLU_valid_LU amesos_klu_l_valid_LU
00108 #define KLU_sort amesos_klu_l_sort
00109 #define KLU_rgrowth amesos_klu_l_rgrowth
00110 #define KLU_rcond amesos_klu_l_rcond
00111 #define KLU_extract amesos_klu_l_extract
00112 #define KLU_condest amesos_klu_l_condest
00113 #define KLU_flops amesos_klu_l_flops
00114 
00115 #else
00116 
00117 #define KLU_scale amesos_klu_scale
00118 #define KLU_solve amesos_klu_solve
00119 #define KLU_tsolve amesos_klu_tsolve
00120 #define KLU_free_numeric amesos_klu_free_numeric
00121 #define KLU_factor amesos_klu_factor
00122 #define KLU_refactor amesos_klu_refactor
00123 #define KLU_kernel_factor amesos_klu_kernel_factor 
00124 #define KLU_lsolve amesos_klu_lsolve
00125 #define KLU_ltsolve amesos_klu_ltsolve
00126 #define KLU_usolve amesos_klu_usolve
00127 #define KLU_utsolve amesos_klu_utsolve
00128 #define KLU_kernel amesos_klu_kernel
00129 #define KLU_valid amesos_klu_valid
00130 #define KLU_valid_LU amesos_klu_valid_LU
00131 #define KLU_sort amesos_klu_sort
00132 #define KLU_rgrowth amesos_klu_rgrowth
00133 #define KLU_rcond amesos_klu_rcond
00134 #define KLU_extract amesos_klu_extract
00135 #define KLU_condest amesos_klu_condest
00136 #define KLU_flops amesos_klu_flops
00137 
00138 #endif
00139 
00140 #endif
00141 
00142 
00143 #ifdef DLONG
00144 
00145 #define KLU_analyze amesos_klu_l_analyze
00146 #define KLU_analyze_given amesos_klu_l_analyze_given
00147 #define KLU_alloc_symbolic amesos_klu_l_alloc_symbolic
00148 #define KLU_free_symbolic amesos_klu_l_free_symbolic
00149 #define KLU_defaults amesos_klu_l_defaults
00150 #define KLU_free amesos_klu_l_free
00151 #define KLU_malloc amesos_klu_l_malloc
00152 #define KLU_realloc amesos_klu_l_realloc
00153 #define KLU_add_size_t amesos_klu_l_add_size_t
00154 #define KLU_mult_size_t amesos_klu_l_mult_size_t
00155 
00156 #define KLU_symbolic klu_l_symbolic
00157 #define KLU_numeric klu_l_numeric
00158 #define KLU_common klu_l_common
00159 
00160 #define BTF_order amesos_btf_l_order
00161 #define BTF_strongcomp amesos_btf_l_strongcomp
00162 
00163 #define AMD_order amesos_amd_l_order
00164 #define COLAMD amesos_colamd_l
00165 #define COLAMD_recommended amesos_colamd_l_recommended
00166 
00167 #else
00168 
00169 #define KLU_analyze amesos_klu_analyze
00170 #define KLU_analyze_given amesos_klu_analyze_given
00171 #define KLU_alloc_symbolic amesos_klu_alloc_symbolic
00172 #define KLU_free_symbolic amesos_klu_free_symbolic
00173 #define KLU_defaults amesos_klu_defaults
00174 #define KLU_free amesos_klu_free
00175 #define KLU_malloc amesos_klu_malloc
00176 #define KLU_realloc amesos_klu_realloc
00177 #define KLU_add_size_t amesos_klu_add_size_t
00178 #define KLU_mult_size_t amesos_klu_mult_size_t
00179 
00180 #define KLU_symbolic klu_symbolic
00181 #define KLU_numeric klu_numeric
00182 #define KLU_common klu_common
00183 
00184 #define BTF_order amesos_btf_order
00185 #define BTF_strongcomp amesos_btf_strongcomp
00186 
00187 #define AMD_order amesos_amd_order
00188 #define COLAMD amesos_colamd
00189 #define COLAMD_recommended amesos_colamd_recommended
00190 
00191 #endif
00192 
00193 
00194 /* -------------------------------------------------------------------------- */
00195 /* Numerical relop macros for correctly handling the NaN case */
00196 /* -------------------------------------------------------------------------- */
00197 
00198 /*
00199 SCALAR_IS_NAN(x):
00200     True if x is NaN.  False otherwise.  The commonly-existing isnan(x)
00201     function could be used, but it's not in Kernighan & Ritchie 2nd edition
00202     (ANSI C).  It may appear in <math.h>, but I'm not certain about
00203     portability.  The expression x != x is true if and only if x is NaN,
00204     according to the IEEE 754 floating-point standard.
00205 
00206 SCALAR_IS_ZERO(x):
00207     True if x is zero.  False if x is nonzero, NaN, or +/- Inf.
00208     This is (x == 0) if the compiler is IEEE 754 compliant.
00209 
00210 SCALAR_IS_NONZERO(x):
00211     True if x is nonzero, NaN, or +/- Inf.  False if x zero.
00212     This is (x != 0) if the compiler is IEEE 754 compliant.
00213 
00214 SCALAR_IS_LTZERO(x):
00215     True if x is < zero or -Inf.  False if x is >= 0, NaN, or +Inf.
00216     This is (x < 0) if the compiler is IEEE 754 compliant.
00217 */
00218 
00219 /* These all work properly, according to the IEEE 754 standard ... except on */
00220 /* a PC with windows.  Works fine in Linux on the same PC... */
00221 #define SCALAR_IS_NAN(x)  ((x) != (x))
00222 #define SCALAR_IS_ZERO(x) ((x) == 0.)
00223 #define SCALAR_IS_NONZERO(x)  ((x) != 0.)
00224 #define SCALAR_IS_LTZERO(x) ((x) < 0.)
00225 
00226 
00227 /* scalar absolute value macro. If x is NaN, the result is NaN: */
00228 #define SCALAR_ABS(x) ((SCALAR_IS_LTZERO (x)) ? -(x) : (x))
00229 
00230 /* print a scalar (avoid printing "-0" for negative zero).  */
00231 #ifdef NPRINT
00232 #define PRINT_SCALAR(a)
00233 #else
00234 #define PRINT_SCALAR(a) \
00235 { \
00236     if (SCALAR_IS_NONZERO (a)) \
00237     { \
00238   PRINTF ((" (%g)", (a))) ; \
00239     } \
00240     else \
00241     { \
00242   PRINTF ((" (0)")) ; \
00243     } \
00244 }
00245 #endif
00246 
00247 /* -------------------------------------------------------------------------- */
00248 /* Real floating-point arithmetic */
00249 /* -------------------------------------------------------------------------- */
00250 
00251 #ifndef COMPLEX
00252 
00253 typedef double Unit ;
00254 #define Entry double
00255 
00256 #define SPLIT(s)            (1)
00257 #define REAL(c)         (c)
00258 #define IMAG(c)         (0.)
00259 #define ASSIGN(c,s1,s2,p,split)     { (c) = (s1)[p] ; }
00260 #define CLEAR(c)        { (c) = 0. ; }
00261 #define CLEAR_AND_INCREMENT(p)      { *p++ = 0. ; }
00262 #define IS_NAN(a)       SCALAR_IS_NAN (a)
00263 #define IS_ZERO(a)        SCALAR_IS_ZERO (a)
00264 #define IS_NONZERO(a)       SCALAR_IS_NONZERO (a)
00265 #define SCALE_DIV(c,s)        { (c) /= (s) ; }
00266 #define SCALE_DIV_ASSIGN(a,c,s)     { a = c / s ; }
00267 #define SCALE(c,s)        { (c) *= (s) ; }
00268 #define ASSEMBLE(c,a)       { (c) += (a) ; }
00269 #define ASSEMBLE_AND_INCREMENT(c,p) { (c) += *p++ ; }
00270 #define DECREMENT(c,a)        { (c) -= (a) ; }
00271 #define MULT(c,a,b)       { (c) = (a) * (b) ; }
00272 #define MULT_CONJ(c,a,b)      { (c) = (a) * (b) ; }
00273 #define MULT_SUB(c,a,b)       { (c) -= (a) * (b) ; }
00274 #define MULT_SUB_CONJ(c,a,b)      { (c) -= (a) * (b) ; }
00275 #define DIV(c,a,b)        { (c) = (a) / (b) ; }
00276 #define RECIPROCAL(c)       { (c) = 1.0 / (c) ; }
00277 #define DIV_CONJ(c,a,b)       { (c) = (a) / (b) ; }
00278 #define APPROX_ABS(s,a)       { (s) = SCALAR_ABS (a) ; }
00279 #define ABS(s,a)        { (s) = SCALAR_ABS (a) ; }
00280 #define PRINT_ENTRY(a)        PRINT_SCALAR (a)
00281 #define CONJ(a,x)       a = x
00282 
00283 /* for flop counts */
00284 #define MULTSUB_FLOPS 2.  /* c -= a*b */
00285 #define DIV_FLOPS 1.  /* c = a/b */
00286 #define ABS_FLOPS 0.  /* c = abs (a) */
00287 #define ASSEMBLE_FLOPS  1.  /* c += a */
00288 #define DECREMENT_FLOPS 1.  /* c -= a */
00289 #define MULT_FLOPS  1.  /* c = a*b */
00290 #define SCALE_FLOPS 1.  /* c = a/s */
00291 
00292 #else
00293 
00294 /* -------------------------------------------------------------------------- */
00295 /* Complex floating-point arithmetic */
00296 /* -------------------------------------------------------------------------- */
00297 
00298 /*
00299     Note:  An alternative to this Double_Complex type would be to use a
00300     struct { double r ; double i ; }.  The problem with that method
00301     (used by the Sun Performance Library, for example) is that ANSI C provides
00302     no guarantee about the layout of a struct.  It is possible that the sizeof
00303     the struct above would be greater than 2 * sizeof (double).  This would
00304     mean that the complex BLAS could not be used.  The method used here avoids
00305     that possibility.  ANSI C *does* guarantee that an array of structs has
00306     the same size as n times the size of one struct.
00307 
00308     The ANSI C99 version of the C language includes a "double _Complex" type.
00309     It should be possible in that case to do the following:
00310 
00311     #define Entry double _Complex
00312 
00313     and remove the Double_Complex struct.  The macros, below, could then be
00314     replaced with instrinsic operators.  Note that the #define Real and
00315     #define Imag should also be removed (they only appear in this file).
00316 
00317     For the MULT, MULT_SUB, MULT_SUB_CONJ, and MULT_CONJ macros,
00318     the output argument c cannot be the same as any input argument.
00319 
00320 */
00321 
00322 typedef struct
00323 {
00324     double component [2] ;  /* real and imaginary parts */
00325 
00326 } Double_Complex ;
00327 
00328 typedef Double_Complex Unit ;
00329 #define Entry Double_Complex
00330 #define Real component [0]
00331 #define Imag component [1]
00332 
00333 /* for flop counts */
00334 #define MULTSUB_FLOPS 8.  /* c -= a*b */
00335 #define DIV_FLOPS 9.  /* c = a/b */
00336 #define ABS_FLOPS 6.  /* c = abs (a), count sqrt as one flop */
00337 #define ASSEMBLE_FLOPS  2.  /* c += a */
00338 #define DECREMENT_FLOPS 2.  /* c -= a */
00339 #define MULT_FLOPS  6.  /* c = a*b */
00340 #define SCALE_FLOPS 2.  /* c = a/s or c = a*s */
00341 
00342 /* -------------------------------------------------------------------------- */
00343 
00344 /* real part of c */
00345 #define REAL(c) ((c).Real)
00346 
00347 /* -------------------------------------------------------------------------- */
00348 
00349 /* imag part of c */
00350 #define IMAG(c) ((c).Imag)
00351 
00352 /* -------------------------------------------------------------------------- */
00353 
00354 /* Return TRUE if a complex number is in split form, FALSE if in packed form */
00355 #define SPLIT(sz) ((sz) != (double *) NULL)
00356 
00357 /* c = (s1) + (s2)*i, if s2 is null, then X is in "packed" format (compatible
00358  * with Entry and ANSI C99 double _Complex type).  */
00359 /*#define ASSIGN(c,s1,s2,p,split) \
00360 { \
00361     if (split) \
00362     { \
00363         (c).Real = (s1)[p] ; \
00364         (c).Imag = (s2)[p] ; \
00365     }  \
00366     else \
00367     { \
00368   (c) = ((Entry *)(s1))[p] ; \
00369     }  \
00370 }*/
00371 
00372 /* -------------------------------------------------------------------------- */
00373 #define CONJ(a, x) \
00374 { \
00375     a.Real = x.Real ; \
00376     a.Imag = -x.Imag ; \
00377 }
00378 
00379 /* c = 0 */
00380 #define CLEAR(c) \
00381 { \
00382     (c).Real = 0. ; \
00383     (c).Imag = 0. ; \
00384 }
00385 
00386 /* -------------------------------------------------------------------------- */
00387 
00388 /* *p++ = 0 */
00389 #define CLEAR_AND_INCREMENT(p) \
00390 { \
00391     p->Real = 0. ; \
00392     p->Imag = 0. ; \
00393     p++ ; \
00394 }
00395 
00396 /* -------------------------------------------------------------------------- */
00397 
00398 /* True if a == 0 */
00399 #define IS_ZERO(a) \
00400     (SCALAR_IS_ZERO ((a).Real) && SCALAR_IS_ZERO ((a).Imag))
00401 
00402 /* -------------------------------------------------------------------------- */
00403 
00404 /* True if a is NaN */
00405 #define IS_NAN(a) \
00406     (SCALAR_IS_NAN ((a).Real) || SCALAR_IS_NAN ((a).Imag))
00407 
00408 /* -------------------------------------------------------------------------- */
00409 
00410 /* True if a != 0 */
00411 #define IS_NONZERO(a) \
00412     (SCALAR_IS_NONZERO ((a).Real) || SCALAR_IS_NONZERO ((a).Imag))
00413 
00414 /* -------------------------------------------------------------------------- */
00415 
00416 /* a = c/s */
00417 #define SCALE_DIV_ASSIGN(a,c,s) \
00418 { \
00419     a.Real = c.Real / s ; \
00420     a.Imag = c.Imag / s ; \
00421 }
00422 
00423 /* c /= s */
00424 #define SCALE_DIV(c,s) \
00425 { \
00426     (c).Real /= (s) ; \
00427     (c).Imag /= (s) ; \
00428 }
00429 
00430 /* -------------------------------------------------------------------------- */
00431 
00432 /* c *= s */
00433 #define SCALE(c,s) \
00434 { \
00435     (c).Real *= (s) ; \
00436     (c).Imag *= (s) ; \
00437 }
00438 
00439 /* -------------------------------------------------------------------------- */
00440 
00441 /* c += a */
00442 #define ASSEMBLE(c,a) \
00443 { \
00444     (c).Real += (a).Real ; \
00445     (c).Imag += (a).Imag ; \
00446 }
00447 
00448 /* -------------------------------------------------------------------------- */
00449 
00450 /* c += *p++ */
00451 #define ASSEMBLE_AND_INCREMENT(c,p) \
00452 { \
00453     (c).Real += p->Real ; \
00454     (c).Imag += p->Imag ; \
00455     p++ ; \
00456 }
00457 
00458 /* -------------------------------------------------------------------------- */
00459 
00460 /* c -= a */
00461 #define DECREMENT(c,a) \
00462 { \
00463     (c).Real -= (a).Real ; \
00464     (c).Imag -= (a).Imag ; \
00465 }
00466 
00467 /* -------------------------------------------------------------------------- */
00468 
00469 /* c = a*b, assert because c cannot be the same as a or b */
00470 #define MULT(c,a,b) \
00471 { \
00472     ASSERT (&(c) != &(a) && &(c) != &(b)) ; \
00473     (c).Real = (a).Real * (b).Real - (a).Imag * (b).Imag ; \
00474     (c).Imag = (a).Imag * (b).Real + (a).Real * (b).Imag ; \
00475 }
00476 
00477 /* -------------------------------------------------------------------------- */
00478 
00479 /* c = a*conjugate(b), assert because c cannot be the same as a or b */
00480 #define MULT_CONJ(c,a,b) \
00481 { \
00482     ASSERT (&(c) != &(a) && &(c) != &(b)) ; \
00483     (c).Real = (a).Real * (b).Real + (a).Imag * (b).Imag ; \
00484     (c).Imag = (a).Imag * (b).Real - (a).Real * (b).Imag ; \
00485 }
00486 
00487 /* -------------------------------------------------------------------------- */
00488 
00489 /* c -= a*b, assert because c cannot be the same as a or b */
00490 #define MULT_SUB(c,a,b) \
00491 { \
00492     ASSERT (&(c) != &(a) && &(c) != &(b)) ; \
00493     (c).Real -= (a).Real * (b).Real - (a).Imag * (b).Imag ; \
00494     (c).Imag -= (a).Imag * (b).Real + (a).Real * (b).Imag ; \
00495 }
00496 
00497 /* -------------------------------------------------------------------------- */
00498 
00499 /* c -= a*conjugate(b), assert because c cannot be the same as a or b */
00500 #define MULT_SUB_CONJ(c,a,b) \
00501 { \
00502     ASSERT (&(c) != &(a) && &(c) != &(b)) ; \
00503     (c).Real -= (a).Real * (b).Real + (a).Imag * (b).Imag ; \
00504     (c).Imag -= (a).Imag * (b).Real - (a).Real * (b).Imag ; \
00505 }
00506 
00507 /* -------------------------------------------------------------------------- */
00508 
00509 /* c = a/b, be careful to avoid underflow and overflow */
00510 #ifdef MATHWORKS
00511 #define DIV(c,a,b) \
00512 { \
00513     (void) utDivideComplex ((a).Real, (a).Imag, (b).Real, (b).Imag, \
00514   &((c).Real), &((c).Imag)) ; \
00515 }
00516 #else
00517 /* This uses ACM Algo 116, by R. L. Smith, 1962. */
00518 /* c can be the same variable as a or b. */
00519 /* Ignore NaN case for double relop br>=bi. */
00520 #define DIV(c,a,b) \
00521 { \
00522     double r, den, ar, ai, br, bi ; \
00523     br = (b).Real ; \
00524     bi = (b).Imag ; \
00525     ar = (a).Real ; \
00526     ai = (a).Imag ; \
00527     if (SCALAR_ABS (br) >= SCALAR_ABS (bi)) \
00528     { \
00529   r = bi / br ; \
00530   den = br + r * bi ; \
00531   (c).Real = (ar + ai * r) / den ; \
00532   (c).Imag = (ai - ar * r) / den ; \
00533     } \
00534     else \
00535     { \
00536   r = br / bi ; \
00537   den = r * br + bi ; \
00538   (c).Real = (ar * r + ai) / den ; \
00539   (c).Imag = (ai * r - ar) / den ; \
00540     } \
00541 }
00542 #endif
00543 
00544 /* -------------------------------------------------------------------------- */
00545 
00546 /* c = 1/c, be careful to avoid underflow and overflow */
00547 /* Not used if MATHWORKS is defined. */
00548 /* This uses ACM Algo 116, by R. L. Smith, 1962. */
00549 /* Ignore NaN case for double relop cr>=ci. */
00550 #define RECIPROCAL(c) \
00551 { \
00552     double r, den, cr, ci ; \
00553     cr = (c).Real ; \
00554     ci = (c).Imag ; \
00555     if (SCALAR_ABS (cr) >= SCALAR_ABS (ci)) \
00556     { \
00557   r = ci / cr ; \
00558   den = cr + r * ci ; \
00559   (c).Real = 1.0 / den ; \
00560   (c).Imag = - r / den ; \
00561     } \
00562     else \
00563     { \
00564   r = cr / ci ; \
00565   den = r * cr + ci ; \
00566   (c).Real = r / den ; \
00567   (c).Imag = - 1.0 / den ; \
00568     } \
00569 }
00570 
00571 
00572 /* -------------------------------------------------------------------------- */
00573 
00574 /* c = a/conjugate(b), be careful to avoid underflow and overflow */
00575 #ifdef MATHWORKS
00576 #define DIV_CONJ(c,a,b) \
00577 { \
00578     (void) utDivideComplex ((a).Real, (a).Imag, (b).Real, (-(b).Imag), \
00579   &((c).Real), &((c).Imag)) ; \
00580 }
00581 #else
00582 /* This uses ACM Algo 116, by R. L. Smith, 1962. */
00583 /* c can be the same variable as a or b. */
00584 /* Ignore NaN case for double relop br>=bi. */
00585 #define DIV_CONJ(c,a,b) \
00586 { \
00587     double r, den, ar, ai, br, bi ; \
00588     br = (b).Real ; \
00589     bi = (b).Imag ; \
00590     ar = (a).Real ; \
00591     ai = (a).Imag ; \
00592     if (SCALAR_ABS (br) >= SCALAR_ABS (bi)) \
00593     { \
00594   r = (-bi) / br ; \
00595   den = br - r * bi ; \
00596   (c).Real = (ar + ai * r) / den ; \
00597   (c).Imag = (ai - ar * r) / den ; \
00598     } \
00599     else \
00600     { \
00601   r = br / (-bi) ; \
00602   den =  r * br - bi; \
00603   (c).Real = (ar * r + ai) / den ; \
00604   (c).Imag = (ai * r - ar) / den ; \
00605     } \
00606 }
00607 #endif
00608 
00609 /* -------------------------------------------------------------------------- */
00610 
00611 /* approximate absolute value, s = |r|+|i| */
00612 #define APPROX_ABS(s,a) \
00613 { \
00614     (s) = SCALAR_ABS ((a).Real) + SCALAR_ABS ((a).Imag) ; \
00615 }
00616 
00617 /* -------------------------------------------------------------------------- */
00618 
00619 /* exact absolute value, s = sqrt (a.real^2 + amag^2) */
00620 #ifdef MATHWORKS
00621 #define ABS(s,a) \
00622 { \
00623     (s) = utFdlibm_hypot ((a).Real, (a).Imag) ; \
00624 }
00625 #else
00626 /* Ignore NaN case for the double relops ar>=ai and ar+ai==ar. */
00627 #define ABS(s,a) \
00628 { \
00629     double r, ar, ai ; \
00630     ar = SCALAR_ABS ((a).Real) ; \
00631     ai = SCALAR_ABS ((a).Imag) ; \
00632     if (ar >= ai) \
00633     { \
00634   if (ar + ai == ar) \
00635   { \
00636       (s) = ar ; \
00637   } \
00638   else \
00639   { \
00640       r = ai / ar ; \
00641       (s) = ar * sqrt (1.0 + r*r) ; \
00642   } \
00643     } \
00644     else \
00645     { \
00646   if (ai + ar == ai) \
00647   { \
00648       (s) = ai ; \
00649   } \
00650   else \
00651   { \
00652       r = ar / ai ; \
00653       (s) = ai * sqrt (1.0 + r*r) ; \
00654   } \
00655     } \
00656 }
00657 #endif
00658 
00659 /* -------------------------------------------------------------------------- */
00660 
00661 /* print an entry (avoid printing "-0" for negative zero).  */
00662 #ifdef NPRINT
00663 #define PRINT_ENTRY(a)
00664 #else
00665 #define PRINT_ENTRY(a) \
00666 { \
00667     if (SCALAR_IS_NONZERO ((a).Real)) \
00668     { \
00669   PRINTF ((" (%g", (a).Real)) ; \
00670     } \
00671     else \
00672     { \
00673   PRINTF ((" (0")) ; \
00674     } \
00675     if (SCALAR_IS_LTZERO ((a).Imag)) \
00676     { \
00677   PRINTF ((" - %gi)", -(a).Imag)) ; \
00678     } \
00679     else if (SCALAR_IS_ZERO ((a).Imag)) \
00680     { \
00681   PRINTF ((" + 0i)")) ; \
00682     } \
00683     else \
00684     { \
00685   PRINTF ((" + %gi)", (a).Imag)) ; \
00686     } \
00687 }
00688 #endif
00689 
00690 /* -------------------------------------------------------------------------- */
00691 
00692 #endif  /* #ifndef COMPLEX */
00693 
00694 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines