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