← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/dgssv.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/dgssv.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 "dsp_defs.h"
 
11
#include "util.h"
 
12
 
 
13
void
 
14
dgssv(SuperMatrix *A, int *perm_c, int *perm_r, SuperMatrix *L,
 
15
      SuperMatrix *U, SuperMatrix *B, int *info )
 
16
{
 
17
/*
 
18
 * Purpose
 
19
 * =======
 
20
 *
 
21
 * DGSSV solves the system of linear equations A*X=B, using the
 
22
 * LU factorization from DGSTRF. 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 = D_; 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 = D_, 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 = D_, Mtype = TRU.
 
100
 *
 
101
 * B       (input/output) SuperMatrix*
 
102
 *         B has types: Stype = DN, Dtype = D_, 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
    double   diag_pivot_thresh = 1.0;
 
125
    double   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 != D_ || 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 != D_ || B->Mtype != GE )
 
140
        *info = -6;
 
141
    if ( *info != 0 ) {
 
142
        i = -(*info);
 
143
        xerbla_("dgssv", &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
        dCreate_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
    dgstrf(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
        dgstrs (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
}