Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_matrixtools_MatrixTimesMatrix_h
23 : #define __PLUMED_matrixtools_MatrixTimesMatrix_h
24 :
25 : #include "core/ActionWithMatrix.h"
26 : #include "core/ParallelTaskManager.h"
27 :
28 : namespace PLMD {
29 : namespace matrixtools {
30 :
31 : namespace helpers {
32 : template <typename T, typename=void>
33 : constexpr bool isDissimilarities=false;
34 : template <typename T>
35 : constexpr bool isDissimilarities<T,std::void_t<typename T::isDissimilarities>> =true;
36 : }
37 :
38 : template <typename T>
39 162 : struct MatrixTimesMatrixInput {
40 : T funcinput;
41 : bool no_thread_gather{false};
42 : bool gatherForceOnColumns{false};
43 : bool secondMatrixIsSparse{false};
44 : std::vector<std::size_t> locations_v;
45 : std::size_t* locations{nullptr};
46 : RequiredMatrixElements outmat;
47 : void update() {
48 35 : locations=locations_v.data();
49 35 : }
50 : #ifdef __PLUMED_HAS_OPENACC
51 : void toACCDevice() const {
52 : #pragma acc enter data copyin(this[0:1],no_thread_gather,gatherForceOnColumns,secondMatrixIsSparse,locations[0:locations_v.size()])
53 : funcinput.toACCDevice();
54 : outmat.toACCDevice();
55 : }
56 : void removeFromACCDevice() const {
57 : funcinput.removeFromACCDevice();
58 : outmat.removeFromACCDevice();
59 : #pragma acc exit data delete(locations[0:locations_v.size()],secondMatrixIsSparse,gatherForceOnColumns,no_thread_gather,this[0:1])
60 : }
61 : #endif //__PLUMED_HAS_OPENACC
62 : };
63 :
64 : class InputVectors {
65 : public:
66 : std::size_t nelem;
67 : std::size_t vec2start;
68 : View<double> arg1;
69 : View<double> arg2;
70 73211 : InputVectors( std::size_t n, double* b ) : nelem(n), vec2start(nelem), arg1(b,n), arg2(b+n,n) {}
71 : };
72 :
73 : template <class CV, typename myPTM=defaultPTM>
74 : class MatrixTimesMatrix : public ActionWithMatrix {
75 : public:
76 : using input_type = MatrixTimesMatrixInput<CV>;
77 : using mytype = MatrixTimesMatrix<CV, myPTM>;
78 : using PTM = typename myPTM::template PTM<mytype>;
79 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
80 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
81 : constexpr static bool isDissimilarities=helpers::isDissimilarities<CV>;
82 : private:
83 : PTM taskmanager;
84 : void getBasicMatrixBookeeping( ArgumentsBookkeeping& argumentsMap ) const ;
85 : public:
86 : static void registerKeywords( Keywords& keys );
87 : explicit MatrixTimesMatrix(const ActionOptions&);
88 : void prepare() override ;
89 : unsigned getNumberOfDerivatives() override;
90 : void calculate() override ;
91 : void getInputData( std::vector<double>& inputdata, ArgumentsBookkeeping& argumentsMap ) const override ;
92 : void getInputData( std::vector<float>& inputdata, ArgumentsBookkeeping& argumentsMap ) const override ;
93 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
94 : static void performTask( std::size_t task_index,
95 : const input_type& actiondata,
96 : ParallelActionsInput& input,
97 : ParallelActionsOutput& output );
98 : static int getNumberOfValuesPerTask( std::size_t task_index,
99 : const input_type& actiondata );
100 : static void getForceIndices( std::size_t task_index,
101 : std::size_t colno,
102 : std::size_t ntotal_force,
103 : const input_type& actiondata,
104 : const ParallelActionsInput& input,
105 : ForceIndexHolder force_indices );
106 : };
107 :
108 : template <class CV, typename myPTM>
109 338 : void MatrixTimesMatrix<CV, myPTM>::registerKeywords( Keywords& keys ) {
110 338 : ActionWithMatrix::registerKeywords(keys);
111 676 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of the output matrix to compute");
112 676 : keys.addInputKeyword("compulsory","ARG","matrix","the label of the two matrices from which the product is calculated");
113 338 : CV::registerKeywords( keys );
114 338 : PTM::registerKeywords( keys );
115 338 : }
116 :
117 : template <class CV, typename myPTM>
118 87 : MatrixTimesMatrix<CV, myPTM>::MatrixTimesMatrix(const ActionOptions&ao):
119 : Action(ao),
120 : ActionWithMatrix(ao),
121 87 : taskmanager(this) {
122 : int nm=getNumberOfMasks();
123 87 : if( nm<0 ) {
124 : nm = 0;
125 : }
126 87 : if( getNumberOfArguments()-nm!=2 ) {
127 0 : error("should be two arguments to this action, a matrix and a vector");
128 : }
129 87 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
130 0 : error("first argument to this action should be a matrix");
131 : }
132 87 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
133 0 : error("second argument to this action should be a matrix");
134 : }
135 87 : if( getPntrToArgument(0)->getShape()[1]!=getPntrToArgument(1)->getShape()[0] ) {
136 0 : error("number of columns in first matrix does not equal number of rows in second matrix");
137 : }
138 87 : std::vector<std::size_t> shape(2);
139 87 : shape[0]=getPntrToArgument(0)->getShape()[0];
140 87 : shape[1]=getPntrToArgument(1)->getShape()[1];
141 87 : addValue( shape );
142 87 : setNotPeriodic();
143 : // With a MASK the output has the sparsity pattern of the mask (set at every
144 : // step by updateBookeepingArrays): do not allocate a dense matrix for it
145 87 : if( nm>0 ) {
146 16 : getPntrToComponent(0)->reshapeMatrixStore( getPntrToArgument(getNumberOfArguments()-1)->getNumberOfColumns() );
147 : } else {
148 71 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
149 : }
150 75 : if( !isDissimilarities && getPntrToArgument(0)->isDerivativeZeroWhenValueIsZero() && getPntrToArgument(1)->isDerivativeZeroWhenValueIsZero() ) {
151 36 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
152 : }
153 :
154 87 : if( nm>0 ) {
155 16 : unsigned iarg = getNumberOfArguments()-1;
156 16 : if( getPntrToArgument(iarg)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
157 0 : error("argument passed to MASK keyword should be a matrix");
158 : }
159 16 : if( getPntrToArgument(iarg)->getShape()[0]!=shape[0] || getPntrToArgument(iarg)->getShape()[1]!=shape[1] ) {
160 0 : error("argument passed to MASK keyword has the wrong shape");
161 : }
162 : }
163 87 : input_type actdata;
164 12 : actdata.funcinput.setup( this, getPntrToArgument(0) );
165 12 : actdata.gatherForceOnColumns = false;
166 87 : actdata.no_thread_gather = no_thread_gather;
167 87 : taskmanager.setActionInput( actdata );
168 87 : }
169 :
170 : template <class CV, typename myPTM>
171 88 : unsigned MatrixTimesMatrix<CV, myPTM>::getNumberOfDerivatives() {
172 88 : return getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getNumberOfStoredValues();
173 : }
174 :
175 : template <class CV, typename myPTM>
176 897 : void MatrixTimesMatrix<CV, myPTM>::prepare() {
177 897 : ActionWithVector::prepare();
178 897 : Value* myval = getPntrToComponent(0);
179 897 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) {
180 880 : return;
181 : }
182 17 : std::vector<std::size_t> shape(2);
183 17 : shape[0]=getPntrToArgument(0)->getShape()[0];
184 17 : shape[1]=getPntrToArgument(1)->getShape()[1];
185 17 : myval->setShape(shape);
186 17 : if( getNumberOfMasks()>0 ) {
187 0 : myval->reshapeMatrixStore( getPntrToArgument(getNumberOfArguments()-1)->getNumberOfColumns() );
188 : } else {
189 17 : myval->reshapeMatrixStore( shape[1] );
190 : }
191 17 : }
192 :
193 : template <class CV, typename myPTM>
194 893 : void MatrixTimesMatrix<CV, myPTM>::calculate() {
195 893 : if( !getPntrToComponent(0)->isDerivativeZeroWhenValueIsZero() ) {
196 749 : if( getPntrToArgument(0)->getNumberOfColumns()<getPntrToArgument(0)->getShape()[1] ) {
197 0 : if( !doNotCalculateDerivatives() ) {
198 0 : error("cannot calculate derivatives for this action with sparse matrices");
199 : } else if(isDissimilarities ) {
200 0 : error("cannot calculate dissimilarities for sparse matrices");
201 : }
202 : }
203 749 : if( getPntrToArgument(1)->getNumberOfColumns()<getPntrToArgument(1)->getShape()[1] ) {
204 6 : if( !doNotCalculateDerivatives() ) {
205 0 : error("cannot calculate derivatives for this action with sparse matrices");
206 : } else if( isDissimilarities ) {
207 0 : error("cannot calculate dissimilarities for sparse matrices");
208 : }
209 : }
210 : }
211 893 : if( getNumberOfMasks()>0 || diagzero ) {
212 128 : updateBookeepingArrays( taskmanager.getActionInput().outmat );
213 : } else {
214 : // Determine the sparsity pattern of the output matrix from the input matrices
215 : Value* arg2 = getPntrToArgument(1);
216 765 : std::vector<unsigned> column_lengths( arg2->getShape()[1], 0 );
217 23552 : for(unsigned i=0; i<arg2->getShape()[0]; ++i) {
218 : unsigned nr = arg2->getRowLength(i);
219 110680 : for(unsigned j=0; j<nr; ++j) {
220 87893 : unsigned ind = arg2->getRowIndex( i, j );
221 87893 : column_lengths[ind]++;
222 : }
223 : }
224 765 : unsigned maxrow = 0;
225 : Value* arg1 = getPntrToArgument(0);
226 5914 : for(unsigned i=0; i<arg1->getShape()[0]; ++i) {
227 5149 : if( arg1->getRowLength(i)==0 ) {
228 0 : continue;
229 : }
230 : unsigned rowlen = 0;
231 567843 : for(unsigned j=0; j<arg2->getShape()[1]; ++j) {
232 562694 : if( column_lengths[j]>0 ) {
233 562694 : rowlen++;
234 : }
235 : }
236 5149 : if( rowlen>maxrow ) {
237 765 : maxrow = rowlen;
238 : }
239 : }
240 765 : Value* myval = getPntrToComponent(0);
241 765 : myval->reshapeMatrixStore( maxrow );
242 : // the store keeps at least one column (see Value::reshapeMatrixStore)
243 765 : maxrow = myval->getNumberOfColumns();
244 5914 : for(unsigned i=0; i<arg1->getShape()[0]; ++i) {
245 5149 : unsigned rstart = i*(maxrow+1);
246 5149 : if( arg1->getRowLength(i)==0 ) {
247 : myval->setMatrixBookeepingElement( rstart, 0 );
248 : continue;
249 : }
250 : unsigned rowlen=0;
251 567843 : for(unsigned j=0; j<arg2->getShape()[1]; ++j) {
252 562694 : if( column_lengths[j]>0 ) {
253 562694 : myval->setMatrixBookeepingElement( rstart + 1 + rowlen, j );
254 562694 : rowlen++;
255 : }
256 : }
257 : myval->setMatrixBookeepingElement( rstart, rowlen );
258 : }
259 765 : copyMatrixBookeepingFromFirstComponent( taskmanager.getActionInput().outmat );
260 : }
261 893 : taskmanager.getActionInput().secondMatrixIsSparse = getPntrToArgument(1)->getNumberOfColumns()<getPntrToArgument(1)->getShape()[1];
262 893 : if( no_thread_gather ) {
263 19 : taskmanager.setupParallelTaskManager( getPntrToArgument(0)->getNumberOfColumns(), 0 );
264 19 : taskmanager.setWorkspaceSize( 4*getPntrToArgument(0)->getNumberOfColumns() );
265 : } else {
266 874 : taskmanager.setupParallelTaskManager( 2*getPntrToArgument(0)->getNumberOfColumns(),
267 : getPntrToArgument(1)->getNumberOfStoredValues() );
268 874 : taskmanager.setWorkspaceSize( 2*getPntrToArgument(0)->getNumberOfColumns() );
269 : }
270 893 : taskmanager.runAllTasks();
271 893 : }
272 :
273 : template <class CV, typename myPTM>
274 893 : void MatrixTimesMatrix<CV, myPTM>::getBasicMatrixBookeeping( ArgumentsBookkeeping& argumentsMap ) const {
275 893 : argumentsMap.nargs = 2;
276 893 : argumentsMap.ranks.resize(2);
277 893 : argumentsMap.ranks[0] = argumentsMap.ranks[1] = 2;
278 893 : argumentsMap.shapestarts.resize(2);
279 893 : argumentsMap.shapestarts[0] = 0;
280 893 : argumentsMap.shapestarts[1] = 2;
281 893 : argumentsMap.ncols.resize(2);
282 893 : argumentsMap.ncols[0] = getPntrToArgument(0)->getNumberOfColumns();
283 893 : argumentsMap.shapedata.resize(4);
284 2679 : for(unsigned i=0; i<2; ++i) {
285 5358 : for(unsigned j=0; j<2; ++j) {
286 3572 : argumentsMap.shapedata[2*i+j] = getPntrToArgument(i)->getShape()[j];
287 : }
288 : }
289 893 : argumentsMap.argstarts.resize(2);
290 893 : argumentsMap.argstarts[0] = 0;
291 893 : argumentsMap.argstarts[1] = getPntrToArgument(0)->getNumberOfStoredValues();
292 :
293 : // Find the lengths of all the columns
294 : Value* arg0 = getPntrToArgument(0);
295 : Value* arg1 = getPntrToArgument(1);
296 893 : argumentsMap.ncols[1] = arg1->getLengthOfLongestColumn();
297 893 : argumentsMap.bookstarts.resize(2);
298 893 : argumentsMap.bookstarts[0] = 0;
299 893 : argumentsMap.bookstarts[1] = arg0->getShape()[0]*(1+argumentsMap.ncols[0]);
300 893 : argumentsMap.booksizes.resize(2);
301 893 : argumentsMap.booksizes[0] = arg0->getShape()[0]*(1+argumentsMap.ncols[0]);
302 893 : argumentsMap.booksizes[1] = arg1->getShape()[1]*(1+argumentsMap.ncols[1]);
303 893 : argumentsMap.bookeeping.resize( argumentsMap.booksizes[0] + argumentsMap.booksizes[1] );
304 156338 : for(unsigned i=0; i<arg0->getShape()[0]; ++i) {
305 : unsigned nr = arg0->getRowLength(i);
306 155445 : unsigned base = i*(argumentsMap.ncols[0]+1);
307 155445 : argumentsMap.bookeeping[base] = nr;
308 4152210 : for(unsigned j=0; j<nr; ++j) {
309 3996765 : argumentsMap.bookeeping[base+1+j] = arg0->getRowIndex( i, j );
310 : }
311 : }
312 155880 : for(unsigned i=0; i<arg1->getShape()[1]; ++i) {
313 154987 : argumentsMap.bookeeping[argumentsMap.bookstarts[1] + i*(argumentsMap.ncols[1]+1)] = 0;
314 : }
315 893 : }
316 :
317 : template <class CV, typename myPTM>
318 893 : void MatrixTimesMatrix<CV, myPTM>::getInputData( std::vector<double>& inputdata, ArgumentsBookkeeping& argumentsMap ) const {
319 893 : getBasicMatrixBookeeping( argumentsMap );
320 893 : std::size_t total_args = getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getShape()[1]*argumentsMap.ncols[1];
321 893 : if( inputdata.size()!=total_args ) {
322 87 : inputdata.resize( total_args );
323 : }
324 893 : getPntrToArgument(0)->assignValues(View{&inputdata[0], getPntrToArgument(0)->getNumberOfStoredValues()});
325 : // The transpose of the input matrix is stored here to make it easier to retrieve the columns
326 :
327 : Value* arg1 = getPntrToArgument(1);
328 893 : unsigned arg_ncol = arg1->getNumberOfColumns();
329 25057 : for(unsigned i=0; i<arg1->getShape()[0]; ++i) {
330 : unsigned nr = arg1->getRowLength(i);
331 3985611 : for(unsigned j=0; j<nr; ++j) {
332 3961447 : unsigned ind = arg1->getRowIndex( i, j );
333 3961447 : unsigned startpos = argumentsMap.bookstarts[1] + ind*(1+argumentsMap.ncols[1]);
334 3961447 : argumentsMap.bookeeping[startpos + 1 + argumentsMap.bookeeping[startpos] ] = i;
335 3961447 : inputdata[ argumentsMap.argstarts[1] + ind*argumentsMap.ncols[1] + argumentsMap.bookeeping[startpos] ] = arg1->get( i*arg_ncol+j, false );
336 3961447 : argumentsMap.bookeeping[startpos]++;
337 : }
338 : }
339 893 : }
340 :
341 : template <class CV, typename myPTM>
342 0 : void MatrixTimesMatrix<CV, myPTM>::getInputData( std::vector<float>& inputdata, ArgumentsBookkeeping& argumentsMap ) const {
343 0 : getBasicMatrixBookeeping( argumentsMap );
344 0 : std::size_t total_args = getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getShape()[1]*argumentsMap.ncols[1];
345 0 : if( inputdata.size()!=total_args ) {
346 0 : inputdata.resize( total_args );
347 : }
348 : Value* arg0 = getPntrToArgument(0);
349 0 : for(unsigned i=0; i<arg0->getNumberOfStoredValues(); ++i) {
350 0 : inputdata[i] = arg0->get(i,false);
351 : }
352 : Value* arg1 = getPntrToArgument(1);
353 0 : unsigned arg_ncol = arg1->getNumberOfColumns();
354 0 : for(unsigned i=0; i<arg1->getShape()[0]; ++i) {
355 : unsigned nr = arg1->getRowLength(i);
356 0 : for(unsigned j=0; j<nr; ++j) {
357 0 : unsigned ind = arg1->getRowIndex( i, j );
358 0 : unsigned startpos = argumentsMap.bookstarts[1] + ind*(1+argumentsMap.ncols[1]);
359 0 : argumentsMap.bookeeping[startpos + 1 + argumentsMap.bookeeping[startpos] ] = i;
360 0 : inputdata[ argumentsMap.argstarts[1] + ind*argumentsMap.ncols[1] + argumentsMap.bookeeping[startpos] ] = arg1->get( i*arg_ncol+j, false );
361 0 : argumentsMap.bookeeping[startpos]++;
362 : }
363 : }
364 0 : }
365 :
366 : template <class CV, typename myPTM>
367 73211 : void MatrixTimesMatrix<CV, myPTM>::performTask( std::size_t task_index,
368 : const input_type& actiondata,
369 : ParallelActionsInput& input,
370 : ParallelActionsOutput& output ) {
371 : // Compute one output row-block for a matrix-times-matrix task.
372 : // Layout notes:
373 : // - arg0/arg1 bookkeeping stores row occupancy and compact index lists.
374 : // - task_index selects one logical row from the first matrix and one row in the output map.
375 : // - output row metadata is encoded as [row_length, col_0, col_1, ...].
376 : // The implementation has two paths:
377 : // (1) sparse second matrix: iterate only explicitly stored nonzeros;
378 : // (2) dense second matrix: iterate contiguous column values.
379 73211 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
380 73211 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
381 : // Number of nonzero/active multiplications for this selected row of arg0.
382 73211 : std::size_t arg0_start = task_index*arg0.ncols;
383 73211 : std::size_t arg0_bookstart = task_index*(1+arg0.ncols);
384 73211 : std::size_t nmult = arg0.bookeeping[arg0_bookstart];
385 : // Output bookkeeping entry for this task row.
386 73211 : std::size_t outrowstart = task_index*(1+actiondata.outmat.ncols);
387 : std::size_t output_rowlen = actiondata.outmat[outrowstart];
388 : // Scratch view used to accumulate the vectorized multiply contributions.
389 : InputVectors vectors( nmult, output.buffer.data() );
390 73211 : if( actiondata.secondMatrixIsSparse ) {
391 605 : if( actiondata.no_thread_gather ) {
392 357 : std::size_t arg1_start = arg1.start + task_index*arg1.ncols;
393 357 : std::size_t arg1_bookstart = task_index*(1+arg1.ncols);
394 2856 : for(unsigned i=0; i<output_rowlen; ++i) {
395 2499 : if( actiondata.gatherForceOnColumns ) {
396 833 : arg0_bookstart = actiondata.outmat[outrowstart+1+i]*(1+arg0.ncols);
397 833 : arg0_start = actiondata.outmat[outrowstart+1+i]*arg0.ncols;
398 : } else {
399 1666 : arg1_bookstart = actiondata.outmat[outrowstart+1+i]*(1+arg1.ncols);
400 1666 : arg1_start = arg1.start + actiondata.outmat[outrowstart+1+i]*arg1.ncols;
401 : }
402 : std::size_t nm = 0;
403 2499 : std::size_t arg1_nelements = arg1.bookeeping[arg1_bookstart];
404 17493 : for(unsigned j=0; j<nmult; ++j) {
405 14994 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
406 59976 : for(unsigned k=0; k<arg1_nelements; ++k) {
407 57834 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
408 12852 : nm++;
409 12852 : break;
410 : }
411 : }
412 : }
413 2499 : if( nm==0 ) {
414 0 : continue ;
415 : }
416 2499 : vectors.nelem = nm;
417 : nm = 0;
418 17493 : for(unsigned j=0; j<nmult; ++j) {
419 14994 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
420 59976 : for(unsigned k=0; k<arg1_nelements; ++k) {
421 57834 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
422 12852 : vectors.arg1[nm] = input.inputdata[ arg0_start + j ];
423 12852 : vectors.arg2[nm] = input.inputdata[ arg1_start + k ];
424 12852 : nm++;
425 12852 : break;
426 : }
427 : }
428 : }
429 2499 : MatrixElementOutput elem( 1, 2*nmult, output.values.data() + i, output.buffer.data() + 2*nmult );
430 2499 : CV::calculate( input.noderiv, actiondata.funcinput, vectors, elem );
431 2499 : if( actiondata.gatherForceOnColumns ) {
432 5117 : for(unsigned k=0; k<nm; ++k) {
433 4284 : output.derivatives[nmult*i + k] = elem.derivs[0][nmult+k];
434 : }
435 : } else {
436 10234 : for(unsigned k=0; k<nm; ++k) {
437 8568 : output.derivatives[nmult*i + k] = elem.derivs[0][k];
438 : }
439 : }
440 : }
441 : } else {
442 2437 : for(unsigned i=0; i<output_rowlen; ++i) {
443 : std::size_t nm = 0;
444 2189 : std::size_t arg1_bookstart = actiondata.outmat[outrowstart+1+i]*(1+arg1.ncols);
445 2189 : std::size_t arg1_nelements = arg1.bookeeping[arg1_bookstart];
446 262332 : for(unsigned j=0; j<nmult; ++j) {
447 260143 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
448 788655 : for(unsigned k=0; k<arg1_nelements; ++k) {
449 538120 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
450 9608 : nm++;
451 9608 : break;
452 : }
453 : }
454 : }
455 2189 : if( nm==0 ) {
456 4 : continue ;
457 : }
458 2185 : vectors.nelem = nm;
459 : nm = 0;
460 2185 : std::size_t arg1_start = arg1.start + actiondata.outmat[outrowstart+1+i]*arg1.ncols;
461 262324 : for(unsigned j=0; j<nmult; ++j) {
462 260139 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
463 788647 : for(unsigned k=0; k<arg1_nelements; ++k) {
464 538116 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
465 9608 : vectors.arg1[nm] = input.inputdata[ arg0_start + j ];
466 9608 : vectors.arg2[nm] = input.inputdata[ arg1_start + k ];
467 9608 : nm++;
468 9608 : break;
469 : }
470 : }
471 : }
472 2185 : MatrixElementOutput elem( 1, 2*nmult, output.values.data() + i, output.derivatives.data() + 2*nmult*i );
473 2185 : CV::calculate( input.noderiv, actiondata.funcinput, vectors, elem );
474 252716 : for(unsigned ii=vectors.nelem; ii<nmult; ++ii) {
475 250531 : elem.derivs[0][ii] = 0;
476 : }
477 : }
478 : }
479 : } else {
480 72606 : if( actiondata.no_thread_gather ) {
481 198 : std::size_t arg1_start = arg1.start + task_index*arg1.ncols;
482 12306 : for(unsigned i=0; i<output_rowlen; ++i) {
483 12108 : if( actiondata.gatherForceOnColumns ) {
484 4036 : arg0_bookstart = actiondata.outmat[outrowstart+1+i]*(1+arg0.ncols);
485 4036 : arg0_start = actiondata.outmat[outrowstart+1+i]*arg0.ncols;
486 : } else {
487 8072 : arg1_start = arg1.start + actiondata.outmat[outrowstart+1+i]*arg1.ncols;
488 : }
489 326748 : for(unsigned j=0; j<nmult; ++j) {
490 314640 : vectors.arg1[j] = input.inputdata[arg0_start + j];
491 314640 : vectors.arg2[j] = input.inputdata[arg1_start + arg0.bookeeping[arg0_bookstart+1+j]];
492 : }
493 12108 : MatrixElementOutput elem( 1, 2*nmult, output.values.data() + i, output.buffer.data() + 2*nmult );
494 12108 : CV::calculate( input.noderiv, actiondata.funcinput, vectors, elem );
495 12108 : if( actiondata.gatherForceOnColumns ) {
496 108916 : for(unsigned k=0; k<nmult; ++k) {
497 104880 : output.derivatives[nmult*i + k] = elem.derivs[0][nmult+k];
498 : }
499 : } else {
500 217832 : for(unsigned k=0; k<nmult; ++k) {
501 209760 : output.derivatives[nmult*i + k] = elem.derivs[0][k];
502 : }
503 : }
504 : }
505 : } else {
506 1850944 : for(unsigned i=0; i<nmult; ++i) {
507 1778536 : vectors.arg1[i] = input.inputdata[arg0_start + i];
508 : }
509 1532149 : for(unsigned i=0; i<output_rowlen; ++i) {
510 1459741 : std::size_t arg1_start = arg1.start + actiondata.outmat[outrowstart+1+i]*arg1.ncols;
511 24459128 : for(unsigned j=0; j<nmult; ++j) {
512 22999387 : vectors.arg2[j] = input.inputdata[arg1_start + arg0.bookeeping[arg0_bookstart+1+j]];
513 : }
514 1459741 : MatrixElementOutput elem( 1, 2*nmult, output.values.data() + i, output.derivatives.data() + 2*nmult*i );
515 1459741 : CV::calculate( input.noderiv, actiondata.funcinput, vectors, elem );
516 : }
517 : }
518 : }
519 73211 : }
520 :
521 : template <class CV, typename myPTM>
522 751 : void MatrixTimesMatrix<CV, myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
523 : Value* arg1 = getPntrToArgument(1);
524 751 : if( arg1->getNumberOfColumns()<arg1->getShape()[1] ) {
525 35 : std::size_t maxcol = arg1->getLengthOfLongestColumn();
526 35 : std::vector<std::size_t>& loc = taskmanager.getActionInput().locations_v;
527 35 : if( loc.size()!=(1+maxcol)*arg1->getShape()[1] ) {
528 3 : loc.resize( (1+maxcol)*arg1->getShape()[1] );
529 : }
530 275 : for(unsigned i=0; i<arg1->getShape()[1]; ++i) {
531 240 : loc[i*(1+maxcol)] = 0;
532 : }
533 35 : unsigned argncols = arg1->getNumberOfColumns();
534 285 : for(unsigned i=0; i<arg1->getShape()[0]; ++i) {
535 : unsigned nr = arg1->getRowLength(i);
536 1680 : for(unsigned j=0; j<nr; ++j) {
537 1430 : unsigned ind = arg1->getRowIndex( i, j );
538 1430 : unsigned startpos = ind*(1+maxcol);
539 1430 : loc[startpos+1+loc[startpos]] = i*argncols + j;
540 1430 : loc[startpos]++;
541 : }
542 : }
543 35 : taskmanager.getActionInput().update();
544 : }
545 751 : taskmanager.applyForces( outforces );
546 751 : if( no_thread_gather ) {
547 19 : getColumnBookeepingArrays( taskmanager.getActionInput().outmat );
548 19 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = true;
549 19 : taskmanager.setNForceScalars( maxcolsize );
550 19 : taskmanager.applyForces( outforces, false );
551 19 : taskmanager.setNForceScalars( getPntrToComponent(0)->getNumberOfColumns() );
552 19 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = false;
553 : }
554 751 : }
555 :
556 : template <class CV, typename myPTM>
557 4081 : int MatrixTimesMatrix<CV, myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
558 : const input_type& actiondata ) {
559 4081 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
560 4081 : return actiondata.outmat[fstart];
561 : }
562 :
563 : template <class CV, typename myPTM>
564 136762 : void MatrixTimesMatrix<CV, myPTM>::getForceIndices( std::size_t task_index,
565 : std::size_t colno,
566 : std::size_t ntotal_force,
567 : const input_type& actiondata,
568 : const ParallelActionsInput& input,
569 : ForceIndexHolder force_indices ) {
570 136762 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
571 136762 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
572 136762 : std::size_t outcol = actiondata.outmat[task_index*(1+actiondata.outmat.ncols)+1+colno];
573 136762 : std::size_t arg0_start = task_index*arg0.ncols;
574 136762 : std::size_t arg0_bookstart = task_index*(1+arg0.ncols);
575 136762 : std::size_t arg1_bookstart = outcol*(1+arg1.ncols);
576 136762 : if( actiondata.gatherForceOnColumns ) {
577 4869 : arg0_start = outcol*arg0.ncols;
578 4869 : arg0_bookstart = outcol*(1+arg0.ncols);
579 4869 : arg1_bookstart = task_index*(1+arg1.ncols);
580 : }
581 136762 : std::size_t nmult = arg0.bookeeping[arg0_bookstart];
582 136762 : if( actiondata.secondMatrixIsSparse ) {
583 : std::size_t nm = 0;
584 2502 : std::size_t arg1_nelements = arg1.bookeeping[arg1_bookstart];
585 2502 : if( !actiondata.no_thread_gather ) {
586 5837 : for(unsigned j=0; j<nmult; ++j) {
587 5001 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
588 19996 : for(unsigned k=0; k<arg1_nelements; ++k) {
589 19281 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
590 4286 : force_indices.indices[0][nm] = arg0_start + j;
591 4286 : force_indices.indices[0][nmult+nm] = arg1.start + actiondata.locations[ arg1_bookstart + 1 + k ];
592 4286 : nm++;
593 4286 : break;
594 : }
595 : }
596 : }
597 : // Right-operand derivatives start at nmult, not at the number of
598 : // matched entries. The gap contains zero left-operand derivatives;
599 : // treating it as thread-unsafe would subtract arg1.start from an
600 : // unused index and access outside the thread-local force buffer.
601 836 : force_indices.threadsafe_derivatives_end[0] = nm==0 ? 0 : nmult;
602 836 : if( nm==0 ) {
603 1 : force_indices.tot_indices[0] = 0;
604 : } else {
605 835 : force_indices.tot_indices[0] = nmult + nm;
606 : }
607 1666 : } else if( actiondata.gatherForceOnColumns ) {
608 5831 : for(unsigned j=0; j<nmult; ++j) {
609 4998 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
610 19992 : for(unsigned k=0; k<arg1_nelements; ++k) {
611 19278 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
612 4284 : force_indices.indices[0][nm] = arg1.start + actiondata.locations[ arg1_bookstart + 1 + k ];
613 4284 : nm++;
614 4284 : break;
615 : }
616 : }
617 : }
618 833 : force_indices.threadsafe_derivatives_end[0] = nm;
619 833 : force_indices.tot_indices[0] = nm;
620 : } else {
621 5831 : for(unsigned j=0; j<nmult; ++j) {
622 4998 : std::size_t arg0_ind = arg0.bookeeping[arg0_bookstart+1+j];
623 19992 : for(unsigned k=0; k<arg1_nelements; ++k) {
624 19278 : if( arg1.bookeeping[arg1_bookstart+1+k]==arg0_ind ) {
625 4284 : force_indices.indices[0][nm] = arg0_start + j;
626 4284 : nm++;
627 4284 : break;
628 : }
629 : }
630 : }
631 833 : force_indices.threadsafe_derivatives_end[0] = nm;
632 833 : force_indices.tot_indices[0] = nm;
633 : }
634 : } else {
635 134260 : if( !actiondata.no_thread_gather ) {
636 1043150 : for(unsigned j=0; j<nmult; ++j) {
637 916962 : force_indices.indices[0][j] = task_index*arg0.ncols + j;
638 916962 : force_indices.indices[0][nmult+j] = arg1.start + arg0.bookeeping[arg0_bookstart+1+j]*arg1.shape[1] + outcol;
639 : }
640 126188 : force_indices.threadsafe_derivatives_end[0] = nmult;
641 126188 : force_indices.tot_indices[0] = 2*nmult;
642 8072 : } else if( actiondata.gatherForceOnColumns ) {
643 108916 : for(unsigned j=0; j<nmult; ++j) {
644 104880 : force_indices.indices[0][j] = arg1.start + arg0.bookeeping[arg0_bookstart+1+j]*arg1.shape[1] + task_index;
645 : }
646 4036 : force_indices.threadsafe_derivatives_end[0] = nmult;
647 4036 : force_indices.tot_indices[0] = nmult;
648 : } else {
649 108916 : for(unsigned j=0; j<nmult; ++j) {
650 104880 : force_indices.indices[0][j] = task_index*arg0.ncols + j;
651 : }
652 4036 : force_indices.threadsafe_derivatives_end[0] = nmult;
653 4036 : force_indices.tot_indices[0] = nmult;
654 : }
655 : }
656 136762 : }
657 :
658 : }
659 : }
660 : #endif
|
