← 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/sgssv.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/sgssv.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
 
#include "ssp_defs.h"
11
 
#include "util.h"
12
 
 
13
 
void
14
 
sgssv(SuperMatrix *A, int *perm_c, int *perm_r, SuperMatrix *L,
15
 
      SuperMatrix *U, SuperMatrix *B, int *info )
16
 
{
17
 
/*
18
 
 * Purpose
19
 
 * =======
20
 
 *
21
 
 * SGSSV solves the system of linear equations A*X=B, using the
22
 
 * LU factorization from SGSTRF. It performs the following steps:
23
 
 *
24
 
 *   1. If A is stored column-wise (A->Stype = 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 = 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
 
 * A       (input) SuperMatrix*
59
 
 *         Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
60
 
 *         of linear equations is A->nrow. Currently, the type of A can be:
61
 
 *         Stype = NC or NR; Dtype = S_; Mtype = GE. In the future, more
62
 
 *         general A will be handled.
63
 
 *
64
 
 * perm_c  (input/output) int*
65
 
 *         If A->Stype = NC, column permutation vector of size A->ncol
66
 
 *         which defines the permutation matrix Pc; perm_c[i] = j means 
67
 
 *         column i of A is in position j in A*Pc.
68
 
 *         On exit, perm_c may be overwritten by the product of the input
69
 
 *         perm_c and a permutation that postorders the elimination tree
70
 
 *         of Pc'*A'*A*Pc; perm_c is not changed if the elimination tree
71
 
 *         is already in postorder.
72
 
 *
73
 
 *         If A->Stype = 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
 
 * perm_r  (output) int*
78
 
 *         If A->Stype = NC, row permutation vector of size A->nrow, 
79
 
 *         which defines the permutation matrix Pr, and is determined 
80
 
 *         by partial pivoting.  perm_r[i] = j means row i of A is in 
81
 
 *         position j in Pr*A.
82
 
 *
83
 
 *         If A->Stype = NR, permutation vector of size A->ncol, which
84
 
 *         determines permutation of rows of transpose(A)
85
 
 *         (columns of A) as described above.
86
 
 *
87
 
 * L       (output) SuperMatrix*
88
 
 *         The factor L from the factorization 
89
 
 *             Pr*A*Pc=L*U              (if A->Stype = NC) or
90
 
 *             Pr*transpose(A)*Pc=L*U   (if A->Stype = NR).
91
 
 *         Uses compressed row subscripts storage for supernodes, i.e.,
92
 
 *         L has types: Stype = SC, Dtype = S_, Mtype = TRLU.
93
 
 *         
94
 
 * U       (output) SuperMatrix*
95
 
 *         The factor U from the factorization 
96
 
 *             Pr*A*Pc=L*U              (if A->Stype = NC) or
97
 
 *             Pr*transpose(A)*Pc=L*U   (if A->Stype = NR).
98
 
 *         Uses column-wise storage scheme, i.e., U has types:
99
 
 *         Stype = NC, Dtype = S_, Mtype = TRU.
100
 
 *
101
 
 * B       (input/output) SuperMatrix*
102
 
 *         B has types: Stype = DN, Dtype = S_, Mtype = GE.
103
 
 *         On entry, the right hand side matrix.
104
 
 *         On exit, the solution matrix if info = 0;
105
 
 *
106
 
 * info    (output) int*
107
 
 *         = 0: successful exit
108
 
 *         > 0: if info = i, and i is
109
 
 *             <= A->ncol: U(i,i) is exactly zero. The factorization has
110
 
 *                been completed, but the factor U is exactly singular,
111
 
 *                so the solution could not be computed.
112
 
 *             > A->ncol: number of bytes allocated when memory allocation
113
 
 *                failure occurred, plus A->ncol.
114
 
 *   
115
 
 */
116
 
    double   t1;        /* Temporary time */
117
 
    char     refact[1], trans[1];
118
 
    DNformat *Bstore;
119
 
    SuperMatrix *AA; /* A in NC format used by the factorization routine.*/
120
 
    SuperMatrix AC; /* Matrix postmultiplied by Pc */
121
 
    int      lwork = 0, *etree, i;
122
 
    
123
 
    /* Set default values for some parameters */
124
 
    float   diag_pivot_thresh = 1.0;
125
 
    float   drop_tol = 0;
126
 
    int      panel_size;     /* panel size */
127
 
    int      relax;          /* no of columns in a relaxed snodes */
128
 
    double   *utime;
129
 
    extern SuperLUStat_t SuperLUStat;
130
 
 
131
 
    /* Test the input parameters ... */
132
 
    *info = 0;
133
 
    Bstore = B->Store;
134
 
    if ( A->nrow != A->ncol || A->nrow < 0 ||
135
 
         (A->Stype != NC && A->Stype != NR) ||
136
 
         A->Dtype != S_ || A->Mtype != GE )
137
 
        *info = -1;
138
 
    else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||
139
 
        B->Stype != DN || B->Dtype != S_ || B->Mtype != GE )
140
 
        *info = -6;
141
 
    if ( *info != 0 ) {
142
 
        i = -(*info);
143
 
        xerbla_("sgssv", &i);
144
 
        return;
145
 
    }
146
 
    
147
 
    *refact = 'N';
148
 
    *trans = 'N';
149
 
    panel_size = sp_ienv(1);
150
 
    relax = sp_ienv(2);
151
 
 
152
 
    StatInit(panel_size, relax);
153
 
    utime = SuperLUStat.utime;
154
 
 
155
 
    /* Convert A to NC format when necessary. */
156
 
    if ( A->Stype == NR ) {
157
 
        NRformat *Astore = A->Store;
158
 
        AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
159
 
        sCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 
160
 
                               Astore->nzval, Astore->colind, Astore->rowptr,
161
 
                               NC, A->Dtype, A->Mtype);
162
 
        *trans = 'T';
163
 
    } else if ( A->Stype == NC ) AA = A;
164
 
 
165
 
    etree = intMalloc(A->ncol);
166
 
 
167
 
    t1 = SuperLU_timer_();
168
 
    sp_preorder(refact, AA, perm_c, etree, &AC);
169
 
    utime[ETREE] = SuperLU_timer_() - t1;
170
 
 
171
 
    /*printf("Factor PA = LU ... relax %d\tw %d\tmaxsuper %d\trowblk %d\n", 
172
 
          relax, panel_size, sp_ienv(3), sp_ienv(4));*/
173
 
    t1 = SuperLU_timer_(); 
174
 
    /* Compute the LU factorization of A. */
175
 
    sgstrf(refact, &AC, diag_pivot_thresh, drop_tol, relax, panel_size,
176
 
           etree, NULL, lwork, perm_r, perm_c, L, U, info);
177
 
    utime[FACT] = SuperLU_timer_() - t1;
178
 
 
179
 
    t1 = SuperLU_timer_();
180
 
    if ( *info == 0 ) {
181
 
        /* Solve the system A*X=B, overwriting B with X. */
182
 
        sgstrs (trans, L, U, perm_r, perm_c, B, info);
183
 
    }
184
 
    utime[SOLVE] = SuperLU_timer_() - t1;
185
 
 
186
 
    SUPERLU_FREE (etree);
187
 
    Destroy_CompCol_Permuted(&AC);
188
 
    if ( A->Stype == NR ) {
189
 
        Destroy_SuperMatrix_Store(AA);
190
 
        SUPERLU_FREE(AA);
191
 
    }
192
 
 
193
 
    PrintStat( &SuperLUStat );
194
 
    StatFree();
195
 
 
196
 
}