← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/dcomplex.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/linsolve/SuperLU/SRC/dcomplex.c
 
1
 
 
2
/*
 
3
 * -- SuperLU routine (version 2.0) --
 
4
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
5
 * and Lawrence Berkeley National Lab.
 
6
 * November 15, 1997
 
7
 *
 
8
 */
 
9
/*
 
10
 * This file defines common arithmetic operations for complex type.
 
11
 */
 
12
#include <math.h>
 
13
#include <stdio.h>
 
14
#include "dcomplex.h"
 
15
 
 
16
 
 
17
/* Complex Division c = a/b */
 
18
void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b)
 
19
{
 
20
    double ratio, den;
 
21
    double abr, abi, cr, ci;
 
22
  
 
23
    if( (abr = b->r) < 0.)
 
24
        abr = - abr;
 
25
    if( (abi = b->i) < 0.)
 
26
        abi = - abi;
 
27
    if( abr <= abi ) {
 
28
        if (abi == 0) {
 
29
            fprintf(stderr, "z_div.c: division by zero");
 
30
            exit (-1);
 
31
        }         
 
32
        ratio = b->r / b->i ;
 
33
        den = b->i * (1 + ratio*ratio);
 
34
        cr = (a->r*ratio + a->i) / den;
 
35
        ci = (a->i*ratio - a->r) / den;
 
36
    } else {
 
37
        ratio = b->i / b->r ;
 
38
        den = b->r * (1 + ratio*ratio);
 
39
        cr = (a->r + a->i*ratio) / den;
 
40
        ci = (a->i - a->r*ratio) / den;
 
41
    }
 
42
    c->r = cr;
 
43
    c->i = ci;
 
44
}
 
45
 
 
46
/* Returns sqrt(z.r^2 + z.i^2) */
 
47
double z_abs(doublecomplex *z)
 
48
{
 
49
    double temp;
 
50
    double real = z->r;
 
51
    double imag = z->i;
 
52
 
 
53
    if (real < 0) real = -real;
 
54
    if (imag < 0) imag = -imag;
 
55
    if (imag > real) {
 
56
        temp = real;
 
57
        real = imag;
 
58
        imag = temp;
 
59
    }
 
60
    if ((real+imag) == real) return(real);
 
61
  
 
62
    temp = imag/real;
 
63
    temp = real*sqrt(1.0 + temp*temp);  /*overflow!!*/
 
64
    return (temp);
 
65
}
 
66
 
 
67
 
 
68
/* Approximates the abs */
 
69
/* Returns abs(z.r) + abs(z.i) */
 
70
double z_abs1(doublecomplex *z)
 
71
{
 
72
    double real = z->r;
 
73
    double imag = z->i;
 
74
  
 
75
    if (real < 0) real = -real;
 
76
    if (imag < 0) imag = -imag;
 
77
 
 
78
    return (real + imag);
 
79
}
 
80
 
 
81
/* Return the exponentiation */
 
82
void z_exp(doublecomplex *r, doublecomplex *z)
 
83
{
 
84
    double expx;
 
85
 
 
86
    expx = exp(z->r);
 
87
    r->r = expx * cos(z->i);
 
88
    r->i = expx * sin(z->i);
 
89
}
 
90
 
 
91
/* Return the complex conjugate */
 
92
void d_cnjg(doublecomplex *r, doublecomplex *z)
 
93
{
 
94
    r->r = z->r;
 
95
    r->i = -z->i;
 
96
}
 
97
 
 
98
/* Return the imaginary part */
 
99
double d_imag(doublecomplex *z)
 
100
{
 
101
    return (z->i);
 
102
}
 
103
 
 
104