← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/scomplex.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/scomplex.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 "scomplex.h"
 
15
 
 
16
 
 
17
/* Complex Division c = a/b */
 
18
void c_div(complex *c, complex *a, complex *b)
 
19
{
 
20
    float ratio, den;
 
21
    float 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
 
 
47
/* Returns sqrt(z.r^2 + z.i^2) */
 
48
double c_abs(complex *z)
 
49
{
 
50
    float temp;
 
51
    float real = z->r;
 
52
    float imag = z->i;
 
53
 
 
54
    if (real < 0) real = -real;
 
55
    if (imag < 0) imag = -imag;
 
56
    if (imag > real) {
 
57
        temp = real;
 
58
        real = imag;
 
59
        imag = temp;
 
60
    }
 
61
    if ((real+imag) == real) return(real);
 
62
  
 
63
    temp = imag/real;
 
64
    temp = real*sqrt(1.0 + temp*temp);  /*overflow!!*/
 
65
    return (temp);
 
66
}
 
67
 
 
68
 
 
69
/* Approximates the abs */
 
70
/* Returns abs(z.r) + abs(z.i) */
 
71
double c_abs1(complex *z)
 
72
{
 
73
    float real = z->r;
 
74
    float imag = z->i;
 
75
  
 
76
    if (real < 0) real = -real;
 
77
    if (imag < 0) imag = -imag;
 
78
 
 
79
    return (real + imag);
 
80
}
 
81
 
 
82
/* Return the exponentiation */
 
83
void c_exp(complex *r, complex *z)
 
84
{
 
85
    float expx;
 
86
 
 
87
    expx = exp(z->r);
 
88
    r->r = expx * cos(z->i);
 
89
    r->i = expx * sin(z->i);
 
90
}
 
91
 
 
92
/* Return the complex conjugate */
 
93
void r_cnjg(complex *r, complex *z)
 
94
{
 
95
    r->r = z->r;
 
96
    r->i = -z->i;
 
97
}
 
98
 
 
99
/* Return the imaginary part */
 
100
double r_imag(complex *z)
 
101
{
 
102
    return (z->i);
 
103
}
 
104
 
 
105