PaStiX Handbook  6.2.1
core_cscalo.c
Go to the documentation of this file.
1 /**
2  *
3  * @file core_cscalo.c
4  *
5  * PaStiX kernel routines
6  *
7  * @copyright 2010-2015 Univ. of Tennessee, Univ. of California Berkeley and
8  * Univ. of Colorado Denver. All rights reserved.
9  * @copyright 2012-2021 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
10  * Univ. Bordeaux. All rights reserved.
11  *
12  * @version 6.2.0
13  * @author Mathieu Faverge
14  * @date 2021-01-03
15  * @generated from /builds/solverstack/pastix/kernels/core_zscalo.c, normal z -> c, Tue Apr 12 09:38:36 2022
16  *
17  **/
18 #include "common.h"
19 #include "blend/solver.h"
20 #include "pastix_ccores.h"
21 #include "cblas.h"
22 #include "kernels_trace.h"
23 
24 /**
25  ******************************************************************************
26  *
27  * @brief Scale a matrix by a diagonal out of place
28  *
29  * Perform the operation: B <- op(A) * D, where A is a general matrix, and D a
30  * diagonal matrix.
31  *
32  *******************************************************************************
33  *
34  * @param[in] trans
35  * @arg PastixNoTrans: No transpose, op( A ) = A;
36  * @arg PastixTrans: Transpose, op( A ) = A;
37  * @arg PastixConjTrans: Conjugate Transpose, op( A ) = conjf(A).
38  *
39  * @param[in] M
40  * Number of rows of the matrix B.
41  * Number of rows of the matrix A.
42  *
43  * @param[in] N
44  * Number of columns of the matrix B.
45  * Number of columns of the matrix A.
46  *
47  * @param[in] A
48  * Matrix of size lda-by-N.
49  *
50  * @param[in] lda
51  * Leading dimension of the array A. lda >= max(1,M).
52  *
53  * @param[in] D
54  * Diagonal matrix of size ldd-by-N.
55  *
56  * @param[in] ldd
57  * Leading dimension of the array D. ldd >= 1.
58  *
59  * @param[inout] B
60  * Matrix of size LDB-by-N.
61  *
62  * @param[in] ldb
63  * Leading dimension of the array B. ldb >= max(1,M)
64  *
65  *******************************************************************************
66  *
67  * @retval PASTIX_SUCCESS successful exit
68  * @retval <0 if -i, the i-th argument had an illegal value
69  * @retval 1, not yet implemented
70  *
71  ******************************************************************************/
72 int
74  pastix_int_t M,
75  pastix_int_t N,
76  const pastix_complex32_t *A,
77  pastix_int_t lda,
78  const pastix_complex32_t *D,
79  pastix_int_t ldd,
80  pastix_complex32_t *B,
81  pastix_int_t ldb )
82 {
83  pastix_complex32_t alpha;
84  pastix_int_t i, j;
85 
86 #if !defined(NDEBUG)
87  if ((trans < PastixNoTrans) ||
88  (trans > PastixConjTrans))
89  {
90  return -1;
91  }
92 
93  if (M < 0) {
94  return -2;
95  }
96  if (N < 0) {
97  return -3;
98  }
99  if ( lda < pastix_imax(1,M) )
100  {
101  return -5;
102  }
103  if ( ldd < 1 )
104  {
105  return -7;
106  }
107  if ( ldb < pastix_imax(1,M) ) {
108  return -9;
109  }
110 #endif
111 
112 #if defined(PRECISION_z) || defined(PRECISION_c)
113  if (trans == PastixConjTrans) {
114  for( j=0; j<N; j++, D += ldd ) {
115  alpha = *D;
116  for( i=0; i<M; i++, B++, A++ ) {
117  *B = conjf(*A) * alpha;
118  }
119  A += lda - M;
120  B += ldb - M;
121  }
122  }
123  else
124 #endif
125  {
126  for( j=0; j<N; j++, D += ldd ) {
127  alpha = *D;
128  for( i=0; i<M; i++, B++, A++ ) {
129  *B = (*A) * alpha;
130  }
131  A += lda - M;
132  B += ldb - M;
133  }
134  }
135 
136  (void)trans;
137  return PASTIX_SUCCESS;
138 }
139 
140 /**
141  *******************************************************************************
142  *
143  * @brief Copy the L term with scaling for the two-terms algorithm
144  *
145  * Performs LD = op(L) * D
146  *
147  *******************************************************************************
148  *
149  * @param[in] trans
150  * @arg PastixNoTrans: No transpose, op( L ) = L;
151  * @arg PastixTrans: Transpose, op( L ) = L;
152  * @arg PastixConjTrans: Conjugate Transpose, op( L ) = conjf(L).
153  *
154  * @param[in] cblk
155  * Pointer to the structure representing the panel to factorize in the
156  * cblktab array. Next column blok must be accessible through cblk[1].
157  *
158  * @param[inout] dataL
159  * The pointer to the correct representation of lower part of the data.
160  * - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
161  * - pastix_lr_block if the block is compressed.
162  *
163  * @param[inout] dataLD
164  * The pointer to the correct representation of LD.
165  * - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
166  * - pastix_lr_block if the block is compressed.
167  *
168  *******************************************************************************/
169 void
171  SolverCblk *cblk,
172  void *dataL,
173  void *dataLD )
174 {
175  const SolverBlok *blok, *lblk;
176  pastix_int_t M, N;
177  pastix_lrblock_t *lrL, *lrLD;
178  pastix_fixdbl_t time;
179  pastix_complex32_t *LD;
180 
181  time = kernel_trace_start( PastixKernelSCALOCblk );
182 
183  N = cblk_colnbr( cblk );
184 
185  blok = cblk->fblokptr + 1; /* Firt off-diagonal block */
186  lblk = cblk[1].fblokptr; /* Next diagonal block */
187 
188  /* if there are off-diagonal supernodes in the column */
189  if ( blok < lblk )
190  {
191  const pastix_complex32_t *L;
192  const pastix_complex32_t *D;
193  pastix_int_t ldl, ldd, ldld;
194 
195  if ( cblk->cblktype & CBLK_COMPRESSED ) {
196  lrL = (pastix_lrblock_t *)dataL;
197  lrLD = (pastix_lrblock_t *)dataLD;
198  D = lrL->u;
199  ldd = N+1;
200 
201  lrL++; lrLD++;
202  for(; blok < lblk; blok++, lrL++, lrLD++) {
203  M = blok_rownbr( blok );
204 
205  assert( lrLD->rk == -1 );
206 
207  /* Copy L in LD */
208  lrLD->rk = lrL->rk;
209  lrLD->rkmax = lrL->rkmax;
210 
211  if ( lrL->rk == -1 ) {
212  assert( M == lrL->rkmax );
213 
214  /* Initialize the workspace */
215  memcpy( lrLD->u, lrL->u, lrL->rkmax * N * sizeof(pastix_complex32_t) );
216  lrLD->v = NULL;
217 
218  L = lrL->u;
219  LD = lrLD->u;
220  }
221  else {
222  /*
223  * Initialize the workspace
224  */
225  memcpy( lrLD->u, lrL->u, M * lrL->rk * sizeof(pastix_complex32_t) );
226  lrLD->v = ((pastix_complex32_t *)lrLD->u) + M * lrL->rk;
227  memcpy( lrLD->v, lrL->v, N * lrL->rkmax * sizeof(pastix_complex32_t) );
228 
229  L = lrL->v;
230  LD = lrLD->v;
231  M = lrLD->rkmax;
232  }
233 
234  ldl = M;
235  ldld = M;
236 
237  /* Compute LD = L * D */
238  core_cscalo( trans, M, N,
239  L, ldl, D, ldd,
240  LD, ldld );
241  }
242  }
243  else if ( cblk->cblktype & CBLK_LAYOUT_2D ) {
244  L = D = (pastix_complex32_t *)dataL;
245  LD = (pastix_complex32_t *)dataLD;
246  ldd = N+1;
247 
248  for(; blok < lblk; blok++) {
249  M = blok_rownbr( blok );
250 
251  /* Compute LD = L * D */
252  core_cscalo( trans, M, N,
253  L + blok->coefind, M, D, ldd,
254  LD + blok->coefind, M );
255  }
256  }
257  else {
258  L = D = (pastix_complex32_t *)dataL;
259  LD = (pastix_complex32_t *)dataLD;
260  ldl = cblk->stride;
261  ldd = cblk->stride+1;
262 
263  M = cblk->stride - N;
264  LD = LD + blok->coefind;
265  ldld = cblk->stride;
266 
267  core_cscalo( trans, M, N, L + blok->coefind, ldl, D, ldd, LD, ldld );
268  }
269  }
270 
271  M = cblk->stride - N;
272  kernel_trace_stop( cblk->fblokptr->inlast, PastixKernelSCALOCblk, M, N, 0, (pastix_fixdbl_t)(M*N), time );
273 }
274 
275 /**
276  *******************************************************************************
277  *
278  * @brief Copy the lower terms of the block with scaling for the two-terms
279  * algorithm.
280  *
281  * Performs B = op(A) * D
282  *
283  *******************************************************************************
284  *
285  * @param[in] trans
286  * @arg PastixNoTrans: No transpose, op( A ) = A;
287  * @arg PastixTrans: Transpose, op( A ) = A;
288  * @arg PastixConjTrans: Conjugate Transpose, op( A ) = conjf(A).
289  *
290  * @param[in] cblk
291  * Pointer to the structure representing the panel to factorize in the
292  * cblktab array. Next column blok must be accessible through cblk[1].
293  *
294  * @param[in] blok_m
295  * Index of the off-diagonal block to be solved in the cblk. All blocks
296  * facing the same cblk, in the current column block will be solved.
297  *
298  * @param[in] dataA
299  * The pointer to the correct representation of data of A.
300  * - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
301  * - pastix_lr_block if the block is compressed.
302  *
303  * @param[in] dataD
304  * The pointer to the correct representation of data of D.
305  * - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
306  * - pastix_lr_block if the block is compressed.
307  *
308  * @param[inout] dataB
309  * The pointer to the correct representation of data of B.
310  * - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
311  * - pastix_lr_block if the block is compressed.
312  *
313  *******************************************************************************/
314 void
316  SolverCblk *cblk,
317  pastix_int_t blok_m,
318  const void *dataA,
319  const void *dataD,
320  void *dataB )
321 {
322  const SolverBlok *fblok, *lblok, *blok;
323  pastix_int_t M, N, ldd, offset, cblk_m;
324  const pastix_complex32_t *lA;
325  pastix_lrblock_t *lrD, *lrB, *lrA;
326  pastix_complex32_t *D, *B, *A;
327  pastix_complex32_t *lB;
328 
329  N = cblk_colnbr( cblk );
330  fblok = cblk[0].fblokptr; /* The diagonal block */
331  lblok = cblk[1].fblokptr; /* The diagonal block of the next cblk */
332  ldd = blok_rownbr( fblok ) + 1;
333 
334  assert( blok_rownbr(fblok) == N );
335  assert( cblk->cblktype & CBLK_LAYOUT_2D );
336 
337  blok = fblok + blok_m;
338  offset = blok->coefind;
339  cblk_m = blok->fcblknm;
340 
341  if ( cblk->cblktype & CBLK_COMPRESSED ) {
342  lrA = (pastix_lrblock_t *)dataA;
343  lrD = (pastix_lrblock_t *)dataD;
344  lrB = (pastix_lrblock_t *)dataB;
345  D = lrD->u;
346  for (; (blok < lblok) && (blok->fcblknm == cblk_m); blok++, lrA++, lrB++) {
347  M = blok_rownbr( blok );
348 
349  /* Copy A in B */
350  lrB->rk = lrA->rk;
351  lrB->rkmax = lrA->rkmax;
352 
353  if ( lrB->rk == -1 ) {
354  assert( M == lrA->rkmax );
355  assert( NULL == lrA->v );
356 
357  /* Initialize the workspace */
358  memcpy( lrB->u, lrA->u, lrA->rkmax * N * sizeof(pastix_complex32_t) );
359  lrB->v = NULL;
360 
361  lA = lrA->u;
362  lB = lrB->u;
363  }
364  else {
365  /*
366  * Initialize the workspace
367  */
368  memcpy( lrB->u, lrA->u, M * lrA->rk * sizeof(pastix_complex32_t) );
369  lrB->v = ((pastix_complex32_t *)lrB->u) + M * lrA->rk;
370  memcpy( lrB->v, lrA->v, N * lrA->rkmax * sizeof(pastix_complex32_t) );
371 
372  lA = lrA->v;
373  lB = lrB->v;
374  M = lrA->rkmax;
375  }
376 
377  /* Compute B = op(A) * D */
378  core_cscalo( trans, M, N,
379  lA, M, D, ldd, lB, M );
380  }
381  }
382  else {
383  A = (pastix_complex32_t *)dataA;
384  D = (pastix_complex32_t *)dataD;
385  B = (pastix_complex32_t *)dataB;
386 
387  for (; (blok < lblok) && (blok->fcblknm == cblk_m); blok++) {
388  lA = A + blok->coefind - offset;
389  lB = B + blok->coefind - offset;
390  M = blok_rownbr(blok);
391 
392  /* Compute B = op(A) * D */
393  core_cscalo( trans, M, N,
394  lA, M, D, ldd, lB, M );
395  }
396  }
397 }
solver.h
blok_rownbr
static pastix_int_t blok_rownbr(const SolverBlok *blok)
Compute the number of rows of a block.
Definition: solver.h:313
cblk_colnbr
static pastix_int_t cblk_colnbr(const SolverCblk *cblk)
Compute the number of columns in a column block.
Definition: solver.h:247
solver_cblk_s::fblokptr
SolverBlok * fblokptr
Definition: solver.h:134
pastix_lrblock_s::v
void * v
Definition: pastix_lowrank.h:116
solver_cblk_s::stride
pastix_int_t stride
Definition: solver.h:135
pastix_lrblock_s::u
void * u
Definition: pastix_lowrank.h:115
solver_cblk_s
Solver column block structure.
Definition: solver.h:127
pastix_trans_t
enum pastix_trans_e pastix_trans_t
Transpostion.
solver_blok_s
Solver block structure.
Definition: solver.h:107
pastix_lrblock_s
The block low-rank structure to hold a matrix in low-rank form.
Definition: pastix_lowrank.h:112
cpucblk_cscalo
void cpucblk_cscalo(pastix_trans_t trans, SolverCblk *cblk, void *dataL, void *dataLD)
Copy the L term with scaling for the two-terms algorithm.
Definition: core_cscalo.c:170
PastixNoTrans
@ PastixNoTrans
Definition: api.h:424
PastixConjTrans
@ PastixConjTrans
Definition: api.h:426
pastix_ccores.h
PASTIX_SUCCESS
@ PASTIX_SUCCESS
Definition: api.h:346
solver_blok_s::coefind
pastix_int_t coefind
Definition: solver.h:114
pastix_lrblock_s::rk
int rk
Definition: pastix_lowrank.h:113
solver_cblk_s::cblktype
int8_t cblktype
Definition: solver.h:130
solver_blok_s::inlast
int8_t inlast
Definition: solver.h:117
cpublok_cscalo
void cpublok_cscalo(pastix_trans_t trans, SolverCblk *cblk, pastix_int_t blok_m, const void *dataA, const void *dataD, void *dataB)
Copy the lower terms of the block with scaling for the two-terms algorithm.
Definition: core_cscalo.c:315
pastix_lrblock_s::rkmax
int rkmax
Definition: pastix_lowrank.h:114
solver_blok_s::fcblknm
pastix_int_t fcblknm
Definition: solver.h:110
core_cscalo
int core_cscalo(pastix_trans_t trans, pastix_int_t M, pastix_int_t N, const pastix_complex32_t *A, pastix_int_t lda, const pastix_complex32_t *D, pastix_int_t ldd, pastix_complex32_t *B, pastix_int_t ldb)
Scale a matrix by a diagonal out of place.
Definition: core_cscalo.c:73