← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/cgssv.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/sparse/SuperLU/SRC/cgssv.c
1
 
 
2
 
 
3
 
/*
4
 
 * -- SuperLU routine (version 3.0) --
5
 
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
6
 
 * and Lawrence Berkeley National Lab.
7
 
 * October 15, 2003
8
 
 *
9
 
 */
10
 
#include "csp_defs.h"
11
 
 
12
 
void
13
 
cgssv(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r,
14
 
      SuperMatrix *L, SuperMatrix *U, SuperMatrix *B,
15
 
      SuperLUStat_t *stat, int *info )
16
 
{
17
 
/*
18
 
 * Purpose
19
 
 * =======
20
 
 *
21
 
 * CGSSV solves the system of linear equations A*X=B, using the
22
 
 * LU factorization from CGSTRF. It performs the following steps:
23
 
 *
24
 
 *   1. If A is stored column-wise (A->Stype = SLU_NC):
25
 
 *
26
 
 *      1.1. Permute the columns of A, forming A*Pc, where Pc
27
 
 *           is a permutation matrix. For more details of this step, 
28
 
 *           see sp_preorder.c.
29
 
 *
30
 
 *      1.2. Factor A as Pr*A*Pc=L*U with the permutation Pr determined
31
 
 *           by Gaussian elimination with partial pivoting.
32
 
 *           L is unit lower triangular with offdiagonal entries
33
 
 *           bounded by 1 in magnitude, and U is upper triangular.
34
 
 *
35
 
 *      1.3. Solve the system of equations A*X=B using the factored
36
 
 *           form of A.
37
 
 *
38
 
 *   2. If A is stored row-wise (A->Stype = SLU_NR), apply the
39
 
 *      above algorithm to the transpose of A:
40
 
 *
41
 
 *      2.1. Permute columns of transpose(A) (rows of A),
42
 
 *           forming transpose(A)*Pc, where Pc is a permutation matrix. 
43
 
 *           For more details of this step, see sp_preorder.c.
44
 
 *
45
 
 *      2.2. Factor A as Pr*transpose(A)*Pc=L*U with the permutation Pr
46
 
 *           determined by Gaussian elimination with partial pivoting.
47
 
 *           L is unit lower triangular with offdiagonal entries
48
 
 *           bounded by 1 in magnitude, and U is upper triangular.
49
 
 *
50
 
 *      2.3. Solve the system of equations A*X=B using the factored
51
 
 *           form of A.
52
 
 *
53
 
 *   See supermatrix.h for the definition of 'SuperMatrix' structure.
54
 
 * 
55
 
 * Arguments
56
 
 * =========
57
 
 *
58
 
 * options (input) superlu_options_t*
59
 
 *         The structure defines the input parameters to control
60
 
 *         how the LU decomposition will be performed and how the
61
 
 *         system will be solved.
62
 
 *
63
 
 * A       (input) SuperMatrix*
64
 
 *         Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
65
 
 *         of linear equations is A->nrow. Currently, the type of A can be:
66
 
 *         Stype = SLU_NC or SLU_NR; Dtype = SLU_C; Mtype = SLU_GE.
67
 
 *         In the future, more general A may be handled.
68
 
 *
69
 
 * perm_c  (input/output) int*
70
 
 *         If A->Stype = SLU_NC, column permutation vector of size A->ncol
71
 
 *         which defines the permutation matrix Pc; perm_c[i] = j means 
72
 
 *         column i of A is in position j in A*Pc.
73
 
 *         If A->Stype = SLU_NR, column permutation vector of size A->nrow
74
 
 *         which describes permutation of columns of transpose(A) 
75
 
 *         (rows of A) as described above.
76
 
 * 
77
 
 *         If options->ColPerm = MY_PERMC or options->Fact = SamePattern or
78
 
 *            options->Fact = SamePattern_SameRowPerm, it is an input argument.
79
 
 *            On exit, perm_c may be overwritten by the product of the input
80
 
 *            perm_c and a permutation that postorders the elimination tree
81
 
 *            of Pc'*A'*A*Pc; perm_c is not changed if the elimination tree
82
 
 *            is already in postorder.
83
 
 *         Otherwise, it is an output argument.
84
 
 * 
85
 
 * perm_r  (input/output) int*
86
 
 *         If A->Stype = SLU_NC, row permutation vector of size A->nrow, 
87
 
 *         which defines the permutation matrix Pr, and is determined 
88
 
 *         by partial pivoting.  perm_r[i] = j means row i of A is in 
89
 
 *         position j in Pr*A.
90
 
 *         If A->Stype = SLU_NR, permutation vector of size A->ncol, which
91
 
 *         determines permutation of rows of transpose(A)
92
 
 *         (columns of A) as described above.
93
 
 *
94
 
 *         If options->RowPerm = MY_PERMR or
95
 
 *            options->Fact = SamePattern_SameRowPerm, perm_r is an
96
 
 *            input argument.
97
 
 *         otherwise it is an output argument.
98
 
 *
99
 
 * L       (output) SuperMatrix*
100
 
 *         The factor L from the factorization 
101
 
 *             Pr*A*Pc=L*U              (if A->Stype = SLU_NC) or
102
 
 *             Pr*transpose(A)*Pc=L*U   (if A->Stype = SLU_NR).
103
 
 *         Uses compressed row subscripts storage for supernodes, i.e.,
104
 
 *         L has types: Stype = SLU_SC, Dtype = SLU_C, Mtype = SLU_TRLU.
105
 
 *         
106
 
 * U       (output) SuperMatrix*
107
 
 *         The factor U from the factorization 
108
 
 *             Pr*A*Pc=L*U              (if A->Stype = SLU_NC) or
109
 
 *             Pr*transpose(A)*Pc=L*U   (if A->Stype = SLU_NR).
110
 
 *         Uses column-wise storage scheme, i.e., U has types:
111
 
 *         Stype = SLU_NC, Dtype = SLU_C, Mtype = SLU_TRU.
112
 
 *
113
 
 * B       (input/output) SuperMatrix*
114
 
 *         B has types: Stype = SLU_DN, Dtype = SLU_C, Mtype = SLU_GE.
115
 
 *         On entry, the right hand side matrix.
116
 
 *         On exit, the solution matrix if info = 0;
117
 
 *
118
 
 * stat   (output) SuperLUStat_t*
119
 
 *        Record the statistics on runtime and floating-point operation count.
120
 
 *        See util.h for the definition of 'SuperLUStat_t'.
121
 
 *
122
 
 * info    (output) int*
123
 
 *         = 0: successful exit
124
 
 *         > 0: if info = i, and i is
125
 
 *             <= A->ncol: U(i,i) is exactly zero. The factorization has
126
 
 *                been completed, but the factor U is exactly singular,
127
 
 *                so the solution could not be computed.
128
 
 *             > A->ncol: number of bytes allocated when memory allocation
129
 
 *                failure occurred, plus A->ncol.
130
 
 *   
131
 
 */
132
 
    DNformat *Bstore;
133
 
    SuperMatrix *AA;/* A in SLU_NC format used by the factorization routine.*/
134
 
    SuperMatrix AC; /* Matrix postmultiplied by Pc */
135
 
    int      lwork = 0, *etree, i;
136
 
    
137
 
    /* Set default values for some parameters */
138
 
    float   drop_tol = 0.;
139
 
    int      panel_size;     /* panel size */
140
 
    int      relax;          /* no of columns in a relaxed snodes */
141
 
    int      permc_spec;
142
 
    trans_t  trans = NOTRANS;
143
 
    double   *utime;
144
 
    double   t; /* Temporary time */
145
 
 
146
 
    /* Test the input parameters ... */
147
 
    *info = 0;
148
 
    Bstore = B->Store;
149
 
    if ( options->Fact != DOFACT ) *info = -1;
150
 
    else if ( A->nrow != A->ncol || A->nrow < 0 ||
151
 
         (A->Stype != SLU_NC && A->Stype != SLU_NR) ||
152
 
         A->Dtype != SLU_C || A->Mtype != SLU_GE )
153
 
        *info = -2;
154
 
    else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||
155
 
        B->Stype != SLU_DN || B->Dtype != SLU_C || B->Mtype != SLU_GE )
156
 
        *info = -7;
157
 
    if ( *info != 0 ) {
158
 
        i = -(*info);
159
 
        xerbla_("cgssv", &i);
160
 
        return;
161
 
    }
162
 
 
163
 
    utime = stat->utime;
164
 
 
165
 
    /* Convert A to SLU_NC format when necessary. */
166
 
    if ( A->Stype == SLU_NR ) {
167
 
        NRformat *Astore = A->Store;
168
 
        AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
169
 
        cCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 
170
 
                               Astore->nzval, Astore->colind, Astore->rowptr,
171
 
                               SLU_NC, A->Dtype, A->Mtype);
172
 
        trans = TRANS;
173
 
    } else {
174
 
        if ( A->Stype == SLU_NC ) AA = A;
175
 
    }
176
 
 
177
 
    t = SuperLU_timer_();
178
 
    /*
179
 
     * Get column permutation vector perm_c[], according to permc_spec:
180
 
     *   permc_spec = NATURAL:  natural ordering 
181
 
     *   permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
182
 
     *   permc_spec = MMD_ATA:  minimum degree on structure of A'*A
183
 
     *   permc_spec = COLAMD:   approximate minimum degree column ordering
184
 
     *   permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
185
 
     */
186
 
    permc_spec = options->ColPerm;
187
 
    if ( permc_spec != MY_PERMC && options->Fact == DOFACT )
188
 
      get_perm_c(permc_spec, AA, perm_c);
189
 
    utime[COLPERM] = SuperLU_timer_() - t;
190
 
 
191
 
    etree = intMalloc(A->ncol);
192
 
 
193
 
    t = SuperLU_timer_();
194
 
    sp_preorder(options, AA, perm_c, etree, &AC);
195
 
    utime[ETREE] = SuperLU_timer_() - t;
196
 
 
197
 
    panel_size = sp_ienv(1);
198
 
    relax = sp_ienv(2);
199
 
 
200
 
    /*printf("Factor PA = LU ... relax %d\tw %d\tmaxsuper %d\trowblk %d\n", 
201
 
          relax, panel_size, sp_ienv(3), sp_ienv(4));*/
202
 
    t = SuperLU_timer_(); 
203
 
    /* Compute the LU factorization of A. */
204
 
    cgstrf(options, &AC, drop_tol, relax, panel_size,
205
 
           etree, NULL, lwork, perm_c, perm_r, L, U, stat, info);
206
 
    utime[FACT] = SuperLU_timer_() - t;
207
 
 
208
 
    t = SuperLU_timer_();
209
 
    if ( *info == 0 ) {
210
 
        /* Solve the system A*X=B, overwriting B with X. */
211
 
        cgstrs (trans, L, U, perm_c, perm_r, B, stat, info);
212
 
    }
213
 
    utime[SOLVE] = SuperLU_timer_() - t;
214
 
 
215
 
    SUPERLU_FREE (etree);
216
 
    Destroy_CompCol_Permuted(&AC);
217
 
    if ( A->Stype == SLU_NR ) {
218
 
        Destroy_SuperMatrix_Store(AA);
219
 
        SUPERLU_FREE(AA);
220
 
    }
221
 
 
222
 
}