← Back to branch summary

~ubuntu-branches/ubuntu/karmic/hypre/karmic

« back to all changes in this revision

Viewing changes to src/FEI_mv/SuperLU/SRC/scolumn_dfs.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-03-20 11:40:12 UTC
  • mfrom: (4.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090320114012-132h6ok9w2r6o609
Tags: 2.4.0b-2
Rebuild against new openmpi.

Show diffs side-by-side

added added

removed removed

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