← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/zgscon.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/superlu/zgscon.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:   zgscon.c
12
 
 * History:     Modified from lapack routines ZGECON.
13
 
 */
14
 
#include <math.h>
15
 
#include "util.h"
16
 
#include "zsp_defs.h"
17
 
 
18
 
void
19
 
zgscon(char *norm, SuperMatrix *L, SuperMatrix *U,
20
 
       double anorm, double *rcond, int *info)
21
 
{
22
 
/*
23
 
    Purpose   
24
 
    =======   
25
 
 
26
 
    ZGSCON 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 ZGETRF.   
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
 
            zgstrf(). Use compressed row subscripts storage for supernodes,
48
 
            i.e., L has types: Stype = SC, Dtype = Z_, Mtype = TRLU.
49
 
 
50
 
    U       (input) SuperMatrix*
51
 
            The factor U from the factorization Pr*A*Pc=L*U as computed by
52
 
            zgstrf(). Use column-wise storage scheme, i.e., U has types:
53
 
            Stype = NC, Dtype = Z_, Mtype = TRU.
54
 
            
55
 
    ANORM   (input) double
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) double*
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
 
    double ainvnm;
73
 
    doublecomplex *work;
74
 
    extern int zrscl_(int *, doublecomplex *, doublecomplex *, int *);
75
 
 
76
 
    extern int zlacon_(int *, doublecomplex *, doublecomplex *, double *, 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 != Z_ || L->Mtype != TRLU)
85
 
         *info = -2;
86
 
    else if (U->nrow < 0 || U->nrow != U->ncol ||
87
 
             U->Stype != NC || U->Dtype != Z_ || U->Mtype != TRU) 
88
 
        *info = -3;
89
 
    if (*info != 0) {
90
 
        i = -(*info);
91
 
        xerbla_("zgscon", &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 = doublecomplexCalloc( 3*L->nrow );
103
 
 
104
 
 
105
 
    if ( !work )
106
 
        ABORT("Malloc fails for work arrays in zgscon.");
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
 
        zlacon_(&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_ztrsv("Lower", "No transpose", "Unit", L, U, &work[0], info);
122
 
 
123
 
            /* Multiply by inv(U). */
124
 
            sp_ztrsv("Upper", "No transpose", "Non-unit", L, U, &work[0],info);
125
 
            
126
 
        } else {
127
 
 
128
 
            /* Multiply by inv(U'). */
129
 
            sp_ztrsv("Upper", "Transpose", "Non-unit", L, U, &work[0], info);
130
 
 
131
 
            /* Multiply by inv(L'). */
132
 
            sp_ztrsv("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
 
} /* zgscon */
145