← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/zgscon.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

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