1
2 /*! @file scomplex.c
3 * \brief Common arithmetic for complex type
4 *
5 * <pre>
6 * -- SuperLU routine (version 2.0) --
7 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
8 * and Lawrence Berkeley National Lab.
9 * November 15, 1997
10 *
11 * This file defines common arithmetic operations for complex type.
12 * </pre>
13 */
14
15 #include <math.h>
16 #include <stdlib.h>
17 #include <stdio.h>
18 #include "slu_scomplex.h"
19
20
21 /*! \brief Complex Division c = a/b */
c_div(complex * c,complex * a,complex * b)22 void c_div(complex *c, complex *a, complex *b)
23 {
24 float ratio, den;
25 float abr, abi, cr, ci;
26
27 if( (abr = b->r) < 0.)
28 abr = - abr;
29 if( (abi = b->i) < 0.)
30 abi = - abi;
31 if( abr <= abi ) {
32 if (abi == 0) {
33 fprintf(stderr, "z_div.c: division by zero\n");
34 exit(-1);
35 }
36 ratio = b->r / b->i ;
37 den = b->i * (1 + ratio*ratio);
38 cr = (a->r*ratio + a->i) / den;
39 ci = (a->i*ratio - a->r) / den;
40 } else {
41 ratio = b->i / b->r ;
42 den = b->r * (1 + ratio*ratio);
43 cr = (a->r + a->i*ratio) / den;
44 ci = (a->i - a->r*ratio) / den;
45 }
46 c->r = cr;
47 c->i = ci;
48 }
49
50
51 /*! \brief Returns sqrt(z.r^2 + z.i^2) */
c_abs(complex * z)52 double c_abs(complex *z)
53 {
54 float temp;
55 float real = z->r;
56 float imag = z->i;
57
58 if (real < 0) real = -real;
59 if (imag < 0) imag = -imag;
60 if (imag > real) {
61 temp = real;
62 real = imag;
63 imag = temp;
64 }
65 if ((real+imag) == real) return(real);
66
67 temp = imag/real;
68 temp = real*sqrt(1.0 + temp*temp); /*overflow!!*/
69 return (temp);
70 }
71
72
73 /*! \brief Approximates the abs. Returns abs(z.r) + abs(z.i) */
c_abs1(complex * z)74 double c_abs1(complex *z)
75 {
76 float real = z->r;
77 float imag = z->i;
78
79 if (real < 0) real = -real;
80 if (imag < 0) imag = -imag;
81
82 return (real + imag);
83 }
84
85 /*! \brief Return the exponentiation */
c_exp(complex * r,complex * z)86 void c_exp(complex *r, complex *z)
87 {
88 float expx;
89
90 expx = exp(z->r);
91 r->r = expx * cos(z->i);
92 r->i = expx * sin(z->i);
93 }
94
95 /*! \brief Return the complex conjugate */
r_cnjg(complex * r,complex * z)96 void r_cnjg(complex *r, complex *z)
97 {
98 r->r = z->r;
99 r->i = -z->i;
100 }
101
102 /*! \brief Return the imaginary part */
r_imag(complex * z)103 double r_imag(complex *z)
104 {
105 return (z->i);
106 }
107
108
109 /*! \brief SIGN functions for complex number. Returns z/abs(z) */
c_sgn(complex * z)110 complex c_sgn(complex *z)
111 {
112 register float t = c_abs(z);
113 register complex retval;
114
115 if (t == 0.0) {
116 retval.r = 1.0, retval.i = 0.0;
117 } else {
118 retval.r = z->r / t, retval.i = z->i / t;
119 }
120
121 return retval;
122 }
123
124 /*! \brief Square-root of a complex number. */
c_sqrt(complex * z)125 complex c_sqrt(complex *z)
126 {
127 complex retval;
128 register float cr, ci, real, imag;
129
130 real = z->r;
131 imag = z->i;
132
133 if ( imag == 0.0 ) {
134 retval.r = sqrt(real);
135 retval.i = 0.0;
136 } else {
137 ci = (sqrt(real*real + imag*imag) - real) / 2.0;
138 ci = sqrt(ci);
139 cr = imag / (2.0 * ci);
140 retval.r = cr;
141 retval.i = ci;
142 }
143
144 return retval;
145 }
146
147
148