SparseLU_column_bmod.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
5 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 /*
12 
13  * NOTE: This file is the modified version of xcolumn_bmod.c file in SuperLU
14 
15  * -- SuperLU routine (version 3.0) --
16  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
17  * and Lawrence Berkeley National Lab.
18  * October 15, 2003
19  *
20  * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
21  *
22  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
23  * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
24  *
25  * Permission is hereby granted to use or copy this program for any
26  * purpose, provided the above notices are retained on all copies.
27  * Permission to modify the code and to distribute modified code is
28  * granted, provided the above notices are retained, and a notice that
29  * the code was modified is included with the above copyright notice.
30  */
31 #ifndef SPARSELU_COLUMN_BMOD_H
32 #define SPARSELU_COLUMN_BMOD_H
33 
34 namespace Eigen {
35 
36 namespace internal {
52 template <typename Scalar, typename StorageIndex>
54  BlockIndexVector segrep, BlockIndexVector repfnz, Index fpanelc, GlobalLU_t& glu)
55 {
56  Index jsupno, k, ksub, krep, ksupno;
57  Index lptr, nrow, isub, irow, nextlu, new_next, ufirst;
58  Index fsupc, nsupc, nsupr, luptr, kfnz, no_zeros;
59  /* krep = representative of current k-th supernode
60  * fsupc = first supernodal column
61  * nsupc = number of columns in a supernode
62  * nsupr = number of rows in a supernode
63  * luptr = location of supernodal LU-block in storage
64  * kfnz = first nonz in the k-th supernodal segment
65  * no_zeros = no lf leading zeros in a supernodal U-segment
66  */
67 
68  jsupno = glu.supno(jcol);
69  // For each nonzero supernode segment of U[*,j] in topological order
70  k = nseg - 1;
71  Index d_fsupc; // distance between the first column of the current panel and the
72  // first column of the current snode
73  Index fst_col; // First column within small LU update
74  Index segsize;
75  for (ksub = 0; ksub < nseg; ksub++)
76  {
77  krep = segrep(k); k--;
78  ksupno = glu.supno(krep);
79  if (jsupno != ksupno )
80  {
81  // outside the rectangular supernode
82  fsupc = glu.xsup(ksupno);
83  fst_col = (std::max)(fsupc, fpanelc);
84 
85  // Distance from the current supernode to the current panel;
86  // d_fsupc = 0 if fsupc > fpanelc
87  d_fsupc = fst_col - fsupc;
88 
89  luptr = glu.xlusup(fst_col) + d_fsupc;
90  lptr = glu.xlsub(fsupc) + d_fsupc;
91 
92  kfnz = repfnz(krep);
93  kfnz = (std::max)(kfnz, fpanelc);
94 
95  segsize = krep - kfnz + 1;
96  nsupc = krep - fst_col + 1;
97  nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc);
98  nrow = nsupr - d_fsupc - nsupc;
99  Index lda = glu.xlusup(fst_col+1) - glu.xlusup(fst_col);
100 
101 
102  // Perform a triangular solver and block update,
103  // then scatter the result of sup-col update to dense
104  no_zeros = kfnz - fst_col;
105  if(segsize==1)
106  LU_kernel_bmod<1>::run(segsize, dense, tempv, glu.lusup, luptr, lda, nrow, glu.lsub, lptr, no_zeros);
107  else
108  LU_kernel_bmod<Dynamic>::run(segsize, dense, tempv, glu.lusup, luptr, lda, nrow, glu.lsub, lptr, no_zeros);
109  } // end if jsupno
110  } // end for each segment
111 
112  // Process the supernodal portion of L\U[*,j]
113  nextlu = glu.xlusup(jcol);
114  fsupc = glu.xsup(jsupno);
115 
116  // copy the SPA dense into L\U[*,j]
117  Index mem;
118  new_next = nextlu + glu.xlsub(fsupc + 1) - glu.xlsub(fsupc);
119  Index offset = internal::first_multiple<Index>(new_next, internal::packet_traits<Scalar>::size) - new_next;
120  if(offset)
121  new_next += offset;
122  while (new_next > glu.nzlumax )
123  {
124  mem = memXpand<ScalarVector>(glu.lusup, glu.nzlumax, nextlu, LUSUP, glu.num_expansions);
125  if (mem) return mem;
126  }
127 
128  for (isub = glu.xlsub(fsupc); isub < glu.xlsub(fsupc+1); isub++)
129  {
130  irow = glu.lsub(isub);
131  glu.lusup(nextlu) = dense(irow);
132  dense(irow) = Scalar(0.0);
133  ++nextlu;
134  }
135 
136  if(offset)
137  {
138  glu.lusup.segment(nextlu,offset).setZero();
139  nextlu += offset;
140  }
141  glu.xlusup(jcol + 1) = StorageIndex(nextlu); // close L\U(*,jcol);
142 
143  /* For more updates within the panel (also within the current supernode),
144  * should start from the first column of the panel, or the first column
145  * of the supernode, whichever is bigger. There are two cases:
146  * 1) fsupc < fpanelc, then fst_col <-- fpanelc
147  * 2) fsupc >= fpanelc, then fst_col <-- fsupc
148  */
149  fst_col = (std::max)(fsupc, fpanelc);
150 
151  if (fst_col < jcol)
152  {
153  // Distance between the current supernode and the current panel
154  // d_fsupc = 0 if fsupc >= fpanelc
155  d_fsupc = fst_col - fsupc;
156 
157  lptr = glu.xlsub(fsupc) + d_fsupc;
158  luptr = glu.xlusup(fst_col) + d_fsupc;
159  nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc); // leading dimension
160  nsupc = jcol - fst_col; // excluding jcol
161  nrow = nsupr - d_fsupc - nsupc;
162 
163  // points to the beginning of jcol in snode L\U(jsupno)
164  ufirst = glu.xlusup(jcol) + d_fsupc;
165  Index lda = glu.xlusup(jcol+1) - glu.xlusup(jcol);
166  MappedMatrixBlock A( &(glu.lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
167  VectorBlock<ScalarVector> u(glu.lusup, ufirst, nsupc);
168  u = A.template triangularView<UnitLower>().solve(u);
169 
170  new (&A) MappedMatrixBlock ( &(glu.lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
171  VectorBlock<ScalarVector> l(glu.lusup, ufirst+nsupc, nrow);
172  l.noalias() -= A * u;
173 
174  } // End if fst_col
175  return 0;
176 }
177 
178 } // end namespace internal
179 } // end namespace Eigen
180 
181 #endif // SPARSELU_COLUMN_BMOD_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
Index column_bmod(const Index jcol, const Index nseg, BlockScalarVector dense, ScalarVector &tempv, BlockIndexVector segrep, BlockIndexVector repfnz, Index fpanelc, GlobalLU_t &glu)
Performs numeric block updates (sup-col) in topological order.
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate offset
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
static EIGEN_DONT_INLINE void run(const Index segsize, BlockScalarVector &dense, ScalarVector &tempv, ScalarVector &lusup, Index &luptr, const Index lda, const Index nrow, IndexVector &lsub, const Index lptr, const Index no_zeros)
Expression of a fixed-size or dynamic-size sub-vector.
static const Line3 l(Rot3(), 1, 1)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
* lda
Definition: eigenvalues.cpp:59
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
The matrix class, also used for vectors and row-vectors.
Convenience specialization of Stride to specify only an outer stride See class Map for some examples...
Definition: Stride.h:106


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:58