1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
/*
* NOTE: This file is the modified version of [s,d,c,z]column_dfs.c file in SuperLU
* -- SuperLU routine (version 2.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November 15, 1997
*
* Copyright (c) 1994 by Xerox Corporation. All rights reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* Permission to modify the code and to distribute modified code is
* granted, provided the above notices are retained, and a notice that
* the code was modified is included with the above copyright notice.
*/
#ifndef SPARSELU_COLUMN_DFS_H
#define SPARSELU_COLUMN_DFS_H
template <typename Scalar, typename StorageIndex> class SparseLUImpl;
namespace Eigen {
namespace internal {
template<typename IndexVector, typename ScalarVector>
struct column_dfs_traits : no_assignment_operator
{
typedef typename ScalarVector::Scalar Scalar;
typedef typename IndexVector::Scalar StorageIndex;
column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl)
: m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
{}
bool update_segrep(Index /*krep*/, Index /*jj*/)
{
return true;
}
void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
{
if (nextl >= m_glu.nzlmax)
m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU;
}
enum { ExpandMem = true };
Index m_jcol;
Index& m_jsuper_ref;
typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& m_glu;
SparseLUImpl<Scalar, StorageIndex>& m_luImpl;
};
/**
* \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
*
* A supernode representative is the last column of a supernode.
* The nonzeros in U[*,j] are segments that end at supernodes representatives.
* The routine returns a list of the supernodal representatives
* in topological order of the dfs that generates them.
* The location of the first nonzero in each supernodal segment
* (supernodal entry location) is also returned.
*
* \param m number of rows in the matrix
* \param jcol Current column
* \param perm_r Row permutation
* \param maxsuper Maximum number of column allowed in a supernode
* \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
* \param lsub_col defines the rhs vector to start the dfs
* \param [in,out] segrep Segment representatives - new segments appended
* \param repfnz First nonzero location in each row
* \param xprune
* \param marker marker[i] == jj, if i was visited during dfs of current column jj;
* \param parent
* \param xplore working array
* \param glu global LU data
* \return 0 success
* > 0 number of bytes allocated when run out of space
*
*/
template <typename Scalar, typename StorageIndex>
Index SparseLUImpl<Scalar,StorageIndex>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,
BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune,
IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
{
Index jsuper = glu.supno(jcol);
Index nextl = glu.xlsub(jcol);
VectorBlock<IndexVector> marker2(marker, 2*m, m);
column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
// For each nonzero in A(*,jcol) do dfs
for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
{
Index krow = lsub_col(k);
lsub_col(k) = emptyIdxLU;
Index kmark = marker2(krow);
// krow was visited before, go to the next nonz;
if (kmark == jcol) continue;
dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
xplore, glu, nextl, krow, traits);
} // for each nonzero ...
Index fsupc;
StorageIndex nsuper = glu.supno(jcol);
StorageIndex jcolp1 = StorageIndex(jcol) + 1;
Index jcolm1 = jcol - 1;
// check to see if j belongs in the same supernode as j-1
if ( jcol == 0 )
{ // Do nothing for column 0
nsuper = glu.supno(0) = 0 ;
}
else
{
fsupc = glu.xsup(nsuper);
StorageIndex jptr = glu.xlsub(jcol); // Not yet compressed
StorageIndex jm1ptr = glu.xlsub(jcolm1);
// Use supernodes of type T2 : see SuperLU paper
if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
// Make sure the number of columns in a supernode doesn't
// exceed threshold
if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU;
/* If jcol starts a new supernode, reclaim storage space in
* glu.lsub from previous supernode. Note we only store
* the subscript set of the first and last columns of
* a supernode. (first for num values, last for pruning)
*/
if (jsuper == emptyIdxLU)
{ // starts a new supernode
if ( (fsupc < jcolm1-1) )
{ // >= 3 columns in nsuper
StorageIndex ito = glu.xlsub(fsupc+1);
glu.xlsub(jcolm1) = ito;
StorageIndex istop = ito + jptr - jm1ptr;
xprune(jcolm1) = istop; // initialize xprune(jcol-1)
glu.xlsub(jcol) = istop;
for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
glu.lsub(ito) = glu.lsub(ifrom);
nextl = ito; // = istop + length(jcol)
}
nsuper++;
glu.supno(jcol) = nsuper;
} // if a new supernode
} // end else: jcol > 0
// Tidy up the pointers before exit
glu.xsup(nsuper+1) = jcolp1;
glu.supno(jcolp1) = nsuper;
xprune(jcol) = StorageIndex(nextl); // Initialize upper bound for pruning
glu.xlsub(jcolp1) = StorageIndex(nextl);
return 0;
}
} // end namespace internal
} // end namespace Eigen
#endif
|
