← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/scolumn_dfs.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/linsolve/SuperLU/SRC/scolumn_dfs.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
 
/*
11
 
  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
12
 
 
13
 
  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
14
 
  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
15
 
 
16
 
  Permission is hereby granted to use or copy this program for any
17
 
  purpose, provided the above notices are retained on all copies.
18
 
  Permission to modify the code and to distribute modified code is
19
 
  granted, provided the above notices are retained, and a notice that
20
 
  the code was modified is included with the above copyright notice.
21
 
*/
22
 
 
23
 
#include "ssp_defs.h"
24
 
 
25
 
/* What type of supernodes we want */
26
 
#define T2_SUPER
27
 
 
28
 
int
29
 
scolumn_dfs(
30
 
           const int  m,         /* in - number of rows in the matrix */
31
 
           const int  jcol,      /* in */
32
 
           int        *perm_r,   /* in */
33
 
           int        *nseg,     /* modified - with new segments appended */
34
 
           int        *lsub_col, /* in - defines the RHS vector to start the dfs */
35
 
           int        *segrep,   /* modified - with new segments appended */
36
 
           int        *repfnz,   /* modified */
37
 
           int        *xprune,   /* modified */
38
 
           int        *marker,   /* modified */
39
 
           int        *parent,   /* working array */
40
 
           int        *xplore,   /* working array */
41
 
           GlobalLU_t *Glu       /* modified */
42
 
           )
43
 
{
44
 
/* 
45
 
 * Purpose
46
 
 * =======
47
 
 *   "column_dfs" performs a symbolic factorization on column jcol, and
48
 
 *   decide the supernode boundary.
49
 
 *
50
 
 *   This routine does not use numeric values, but only use the RHS 
51
 
 *   row indices to start the dfs.
52
 
 *
53
 
 *   A supernode representative is the last column of a supernode.
54
 
 *   The nonzeros in U[*,j] are segments that end at supernodal
55
 
 *   representatives. The routine returns a list of such supernodal 
56
 
 *   representatives in topological order of the dfs that generates them.
57
 
 *   The location of the first nonzero in each such supernodal segment
58
 
 *   (supernodal entry location) is also returned.
59
 
 *
60
 
 * Local parameters
61
 
 * ================
62
 
 *   nseg: no of segments in current U[*,j]
63
 
 *   jsuper: jsuper=EMPTY if column j does not belong to the same
64
 
 *      supernode as j-1. Otherwise, jsuper=nsuper.
65
 
 *
66
 
 *   marker2: A-row --> A-row/col (0/1)
67
 
 *   repfnz: SuperA-col --> PA-row
68
 
 *   parent: SuperA-col --> SuperA-col
69
 
 *   xplore: SuperA-col --> index to L-structure
70
 
 *
71
 
 * Return value
72
 
 * ============
73
 
 *     0  success;
74
 
 *   > 0  number of bytes allocated when run out of space.
75
 
 *
76
 
 */
77
 
    int     jcolp1, jcolm1, jsuper, nsuper, nextl;
78
 
    int     k, krep, krow, kmark, kperm;
79
 
    int     *marker2;           /* Used for small panel LU */
80
 
    int     fsupc;              /* First column of a snode */
81
 
    int     myfnz;              /* First nonz column of a U-segment */
82
 
    int     chperm, chmark, chrep, kchild;
83
 
    int     xdfs, maxdfs, kpar, oldrep;
84
 
    int     jptr, jm1ptr;
85
 
    int     ito, ifrom, istop;  /* Used to compress row subscripts */
86
 
    int     mem_error;
87
 
    int     *xsup, *supno, *lsub, *xlsub;
88
 
    int     nzlmax;
89
 
    static  int  first = 1, maxsuper;
90
 
    
91
 
    xsup    = Glu->xsup;
92
 
    supno   = Glu->supno;
93
 
    lsub    = Glu->lsub;
94
 
    xlsub   = Glu->xlsub;
95
 
    nzlmax  = Glu->nzlmax;
96
 
 
97
 
    if ( first ) {
98
 
        maxsuper = sp_ienv(3);
99
 
        first = 0;
100
 
    }
101
 
    jcolp1  = jcol + 1;
102
 
    jcolm1  = jcol - 1;
103
 
    nsuper  = supno[jcol];
104
 
    jsuper  = nsuper;
105
 
    nextl   = xlsub[jcol];
106
 
    marker2 = &marker[2*m];
107
 
 
108
 
 
109
 
    /* For each nonzero in A[*,jcol] do dfs */
110
 
    for (k = 0; lsub_col[k] != EMPTY; k++) {
111
 
 
112
 
        krow = lsub_col[k];
113
 
        lsub_col[k] = EMPTY;
114
 
        kmark = marker2[krow];          
115
 
 
116
 
        /* krow was visited before, go to the next nonz */
117
 
        if ( kmark == jcol ) continue; 
118
 
 
119
 
        /* For each unmarked nbr krow of jcol
120
 
         *      krow is in L: place it in structure of L[*,jcol]
121
 
         */
122
 
        marker2[krow] = jcol;
123
 
        kperm = perm_r[krow];
124
 
 
125
 
        if ( kperm == EMPTY ) {
126
 
            lsub[nextl++] = krow;       /* krow is indexed into A */
127
 
            if ( nextl >= nzlmax ) {
128
 
                if ( mem_error = sLUMemXpand(jcol, nextl, LSUB, &nzlmax, Glu) )
129
 
                    return (mem_error);
130
 
                lsub = Glu->lsub;
131
 
            }
132
 
            if ( kmark != jcolm1 ) jsuper = EMPTY;/* Row index subset testing */
133
 
        } else {
134
 
            /*  krow is in U: if its supernode-rep krep
135
 
             *  has been explored, update repfnz[*]
136
 
             */
137
 
            krep = xsup[supno[kperm]+1] - 1;
138
 
            myfnz = repfnz[krep];
139
 
 
140
 
            if ( myfnz != EMPTY ) {     /* Visited before */
141
 
                if ( myfnz > kperm ) repfnz[krep] = kperm;
142
 
                /* continue; */
143
 
            }
144
 
            else {
145
 
                /* Otherwise, perform dfs starting at krep */
146
 
                oldrep = EMPTY;
147
 
                parent[krep] = oldrep;
148
 
                repfnz[krep] = kperm;
149
 
                xdfs = xlsub[krep];
150
 
                maxdfs = xprune[krep];
151
 
 
152
 
                do {
153
 
                    /* 
154
 
                     * For each unmarked kchild of krep 
155
 
                     */
156
 
                    while ( xdfs < maxdfs ) {
157
 
 
158
 
                        kchild = lsub[xdfs];
159
 
                        xdfs++;
160
 
                        chmark = marker2[kchild];
161
 
 
162
 
                        if ( chmark != jcol ) { /* Not reached yet */
163
 
                            marker2[kchild] = jcol;
164
 
                            chperm = perm_r[kchild];
165
 
 
166
 
                            /* Case kchild is in L: place it in L[*,k] */
167
 
                            if ( chperm == EMPTY ) {
168
 
                                lsub[nextl++] = kchild;
169
 
                                if ( nextl >= nzlmax ) {
170
 
                                    if ( mem_error =
171
 
                                         sLUMemXpand(jcol,nextl,LSUB,&nzlmax,Glu) )
172
 
                                        return (mem_error);
173
 
                                    lsub = Glu->lsub;
174
 
                                }
175
 
                                if ( chmark != jcolm1 ) jsuper = EMPTY;
176
 
                            } else {
177
 
                                /* Case kchild is in U: 
178
 
                                 *   chrep = its supernode-rep. If its rep has 
179
 
                                 *   been explored, update its repfnz[*]
180
 
                                 */
181
 
                                chrep = xsup[supno[chperm]+1] - 1;
182
 
                                myfnz = repfnz[chrep];
183
 
                                if ( myfnz != EMPTY ) { /* Visited before */
184
 
                                    if ( myfnz > chperm )
185
 
                                        repfnz[chrep] = chperm;
186
 
                                } else {
187
 
                                    /* Continue dfs at super-rep of kchild */
188
 
                                    xplore[krep] = xdfs;        
189
 
                                    oldrep = krep;
190
 
                                    krep = chrep; /* Go deeper down G(L^t) */
191
 
                                    parent[krep] = oldrep;
192
 
                                    repfnz[krep] = chperm;
193
 
                                    xdfs = xlsub[krep];     
194
 
                                    maxdfs = xprune[krep];
195
 
                                } /* else */
196
 
 
197
 
                           } /* else */
198
 
 
199
 
                        } /* if */
200
 
 
201
 
                    } /* while */
202
 
 
203
 
                    /* krow has no more unexplored nbrs;
204
 
                     *    place supernode-rep krep in postorder DFS.
205
 
                     *    backtrack dfs to its parent
206
 
                     */
207
 
                    segrep[*nseg] = krep;
208
 
                    ++(*nseg);
209
 
                    kpar = parent[krep]; /* Pop from stack, mimic recursion */
210
 
                    if ( kpar == EMPTY ) break; /* dfs done */
211
 
                    krep = kpar;
212
 
                    xdfs = xplore[krep];
213
 
                    maxdfs = xprune[krep];
214
 
 
215
 
                } while ( kpar != EMPTY );      /* Until empty stack */
216
 
 
217
 
            } /* else */
218
 
 
219
 
        } /* else */
220
 
 
221
 
    } /* for each nonzero ... */
222
 
 
223
 
    /* Check to see if j belongs in the same supernode as j-1 */
224
 
    if ( jcol == 0 ) { /* Do nothing for column 0 */
225
 
        nsuper = supno[0] = 0;
226
 
    } else {
227
 
        fsupc = xsup[nsuper];
228
 
        jptr = xlsub[jcol];     /* Not compressed yet */
229
 
        jm1ptr = xlsub[jcolm1];
230
 
 
231
 
#ifdef T2_SUPER
232
 
        if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = EMPTY;
233
 
#endif
234
 
        /* Make sure the number of columns in a supernode doesn't
235
 
           exceed threshold. */
236
 
        if ( jcol - fsupc >= maxsuper ) jsuper = EMPTY;
237
 
 
238
 
        /* If jcol starts a new supernode, reclaim storage space in
239
 
         * lsub from the previous supernode. Note we only store
240
 
         * the subscript set of the first and last columns of
241
 
         * a supernode. (first for num values, last for pruning)
242
 
         */
243
 
        if ( jsuper == EMPTY ) {        /* starts a new supernode */
244
 
            if ( (fsupc < jcolm1-1) ) { /* >= 3 columns in nsuper */
245
 
#ifdef CHK_COMPRESS
246
 
                printf("  Compress lsub[] at super %d-%d\n", fsupc, jcolm1);
247
 
#endif
248
 
                ito = xlsub[fsupc+1];
249
 
                xlsub[jcolm1] = ito;
250
 
                istop = ito + jptr - jm1ptr;
251
 
                xprune[jcolm1] = istop; /* Initialize xprune[jcol-1] */
252
 
                xlsub[jcol] = istop;
253
 
                for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
254
 
                    lsub[ito] = lsub[ifrom];
255
 
                nextl = ito;            /* = istop + length(jcol) */
256
 
            }
257
 
            nsuper++;
258
 
            supno[jcol] = nsuper;
259
 
        } /* if a new supernode */
260
 
 
261
 
    }   /* else: jcol > 0 */ 
262
 
    
263
 
    /* Tidy up the pointers before exit */
264
 
    xsup[nsuper+1] = jcolp1;
265
 
    supno[jcolp1]  = nsuper;
266
 
    xprune[jcol]   = nextl;     /* Initialize upper bound for pruning */
267
 
    xlsub[jcolp1]  = nextl;
268
 
 
269
 
    return 0;
270
 
}