← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/cgscon.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/sandbox/pysparse/superlu/cgscon.c
 
1
 
 
2
 
 
3
/*
 
4
 * -- SuperLU routine (version 2.0) --
 
5
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
6
 * and Lawrence Berkeley National Lab.
 
7
 * November 15, 1997
 
8
 *
 
9
 */
 
10
/*
 
11
 * File name:   cgscon.c
 
12
 * History:     Modified from lapack routines CGECON.
 
13
 */
 
14
#include <math.h>
 
15
#include "util.h"
 
16
#include "csp_defs.h"
 
17
 
 
18
void
 
19
cgscon(char *norm, SuperMatrix *L, SuperMatrix *U,
 
20
       float anorm, float *rcond, int *info)
 
21
{
 
22
/*
 
23
    Purpose   
 
24
    =======   
 
25
 
 
26
    CGSCON estimates the reciprocal of the condition number of a general 
 
27
    real matrix A, in either the 1-norm or the infinity-norm, using   
 
28
    the LU factorization computed by CGETRF.   
 
29
 
 
30
    An estimate is obtained for norm(inv(A)), and the reciprocal of the   
 
31
    condition number is computed as   
 
32
       RCOND = 1 / ( norm(A) * norm(inv(A)) ).   
 
33
 
 
34
    See supermatrix.h for the definition of 'SuperMatrix' structure.
 
35
 
 
36
    Arguments   
 
37
    =========   
 
38
 
 
39
    NORM    (input) char*
 
40
            Specifies whether the 1-norm condition number or the   
 
41
            infinity-norm condition number is required:   
 
42
            = '1' or 'O':  1-norm;   
 
43
            = 'I':         Infinity-norm.
 
44
            
 
45
    L       (input) SuperMatrix*
 
46
            The factor L from the factorization Pr*A*Pc=L*U as computed by
 
47
            cgstrf(). Use compressed row subscripts storage for supernodes,
 
48
            i.e., L has types: Stype = SC, Dtype = C_, Mtype = TRLU.
 
49
 
 
50
    U       (input) SuperMatrix*
 
51
            The factor U from the factorization Pr*A*Pc=L*U as computed by
 
52
            cgstrf(). Use column-wise storage scheme, i.e., U has types:
 
53
            Stype = NC, Dtype = C_, Mtype = TRU.
 
54
            
 
55
    ANORM   (input) float
 
56
            If NORM = '1' or 'O', the 1-norm of the original matrix A.   
 
57
            If NORM = 'I', the infinity-norm of the original matrix A.
 
58
            
 
59
    RCOND   (output) float*
 
60
            The reciprocal of the condition number of the matrix A,   
 
61
            computed as RCOND = 1/(norm(A) * norm(inv(A))).
 
62
            
 
63
    INFO    (output) int*
 
64
            = 0:  successful exit   
 
65
            < 0:  if INFO = -i, the i-th argument had an illegal value   
 
66
 
 
67
    ===================================================================== 
 
68
*/
 
69
 
 
70
    /* Local variables */
 
71
    int    kase, kase1, onenrm, i;
 
72
    float ainvnm;
 
73
    complex *work;
 
74
    extern int crscl_(int *, complex *, complex *, int *);
 
75
 
 
76
    extern int clacon_(int *, complex *, complex *, float *, int *);
 
77
 
 
78
    
 
79
    /* Test the input parameters. */
 
80
    *info = 0;
 
81
    onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O");
 
82
    if (! onenrm && ! lsame_(norm, "I")) *info = -1;
 
83
    else if (L->nrow < 0 || L->nrow != L->ncol ||
 
84
             L->Stype != SC || L->Dtype != C_ || L->Mtype != TRLU)
 
85
         *info = -2;
 
86
    else if (U->nrow < 0 || U->nrow != U->ncol ||
 
87
             U->Stype != NC || U->Dtype != C_ || U->Mtype != TRU) 
 
88
        *info = -3;
 
89
    if (*info != 0) {
 
90
        i = -(*info);
 
91
        xerbla_("cgscon", &i);
 
92
        return;
 
93
    }
 
94
 
 
95
    /* Quick return if possible */
 
96
    *rcond = 0.;
 
97
    if ( L->nrow == 0 || U->nrow == 0) {
 
98
        *rcond = 1.;
 
99
        return;
 
100
    }
 
101
 
 
102
    work = complexCalloc( 3*L->nrow );
 
103
 
 
104
 
 
105
    if ( !work )
 
106
        ABORT("Malloc fails for work arrays in cgscon.");
 
107
    
 
108
    /* Estimate the norm of inv(A). */
 
109
    ainvnm = 0.;
 
110
    if ( onenrm ) kase1 = 1;
 
111
    else kase1 = 2;
 
112
    kase = 0;
 
113
 
 
114
    do {
 
115
        clacon_(&L->nrow, &work[L->nrow], &work[0], &ainvnm, &kase);
 
116
 
 
117
        if (kase == 0) break;
 
118
 
 
119
        if (kase == kase1) {
 
120
            /* Multiply by inv(L). */
 
121
            sp_ctrsv("Lower", "No transpose", "Unit", L, U, &work[0], info);
 
122
 
 
123
            /* Multiply by inv(U). */
 
124
            sp_ctrsv("Upper", "No transpose", "Non-unit", L, U, &work[0],info);
 
125
            
 
126
        } else {
 
127
 
 
128
            /* Multiply by inv(U'). */
 
129
            sp_ctrsv("Upper", "Transpose", "Non-unit", L, U, &work[0], info);
 
130
 
 
131
            /* Multiply by inv(L'). */
 
132
            sp_ctrsv("Lower", "Transpose", "Unit", L, U, &work[0], info);
 
133
            
 
134
        }
 
135
 
 
136
    } while ( kase != 0 );
 
137
 
 
138
    /* Compute the estimate of the reciprocal condition number. */
 
139
    if (ainvnm != 0.) *rcond = (1. / ainvnm) / anorm;
 
140
 
 
141
    SUPERLU_FREE (work);
 
142
    return;
 
143
 
 
144
} /* cgscon */
 
145