1 /* Siconos is a program dedicated to modeling, simulation and control
2 * of non smooth dynamical systems.
3 *
4 * Copyright 2021 INRIA.
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License");
7 * you may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 */
18
19 #include "CSparseMatrix_internal.h"
20 #include "SparseBlockMatrix.h"
21 #include <assert.h> // for assert
22 #include <math.h> // for fabs, NAN
23 #include <stdio.h> // for size_t, fprintf, printf, fscanf, NULL
24 #include <stdlib.h> // for malloc, free, exit, realloc, calloc
25 #include <stdint.h> // for SIZE_MAX
26 #include <string.h> // for memcpy
27 #include "NumericsArrays.h" // for NA_merge_and_sort_sorted_arrays
28 #include "SiconosBlas.h" // for cblas_dscal, cblas_dgemv, CblasNoTrans, max
29 #include "SiconosCompat.h" // for SN_SIZE_T_Fn
30 #include "SiconosLapack.h" // for lapack_int, DGETRF, DGETRI
31 /* #define DEBUG_NOCOLOR 1 */
32 /* #define DEBUG_STDOUT 1 */
33 /* #define DEBUG_MESSAGES 1 */
34 #include "siconos_debug.h" // for DEBUG_PRINTF, DEBUG_END, DEBUG_BEGIN
35 #include "numerics_verbose.h" // for CHECK_IO, numerics_error, numerics_war...
36 #include "op3x3.h" // for mvp3x3, mvp_alpha3x3
37 #include "NSSTools.h" // for min, max
38
39
40
41 #ifdef DEBUG_MESSAGES
42 #include "NumericsMatrix.h"
43 #endif
44
45 //#define VERBOSE_DEBUG
46
47 /* Clang 3.6 seems to be at odd with the C99 and C11 std on the conversion from
48 * double to _Bool, cf paragraph "6.3.1.2 Boolean type" in the latest C99 or C11 draft */
49 #ifdef __clang__
50 #if (__clang_major__ == 3 && __clang_minor__ == 6)
51 #pragma clang diagnostic ignored "-Wfloat-conversion"
52 #endif
53 #endif
54
SBM_inc_version(SparseBlockStructuredMatrix * M)55 static void SBM_inc_version(SparseBlockStructuredMatrix* M)
56 {
57 NDV_inc(&(M->version));
58 }
59
SBM_null(SparseBlockStructuredMatrix * sbm)60 void SBM_null(SparseBlockStructuredMatrix* sbm)
61 {
62 sbm->nbblocks = 0;
63 sbm->block = NULL;
64 sbm->blocknumber0 = 0;
65 sbm->blocknumber1 = 0;
66 sbm->blocksize0 = NULL;
67 sbm->blocksize1 = NULL;
68 sbm->filled1 = 0;
69 sbm->filled2 = 0;
70 sbm->index1_data = NULL;
71 sbm->index2_data = NULL;
72 sbm->diagonal_blocks = NULL;
73
74 NDV_reset(&(sbm->version));
75 }
76
SBM_new(void)77 SparseBlockStructuredMatrix* SBM_new(void)
78 {
79 SparseBlockStructuredMatrix* sbm = (SparseBlockStructuredMatrix*)
80 malloc(sizeof(SparseBlockStructuredMatrix));
81
82 SBM_null(sbm);
83
84 return sbm;
85 }
86
SBM_clear(SparseBlockStructuredMatrix * sbm)87 void SBM_clear(SparseBlockStructuredMatrix *sbm)
88 {
89 /* Free memory for SparseBlockStructuredMatrix */
90 /* Warning: nothing is done to check if memory has really been
91 allocated for each component or if it was only pointer links. Note
92 that when used from the Kernel, memory is not directly allocated
93 for such structures and this function must not be called. See in
94 kernel/src/simulationsTools/SparseBlockMatrix.cpp for details on
95 the way the structure is filled in.
96 */
97 assert(sbm);
98
99 if(sbm->blocksize0)
100 {
101
102 free(sbm->blocksize0);
103 if(sbm->blocksize0 == sbm->blocksize1)
104 {
105 sbm->blocksize1 = NULL ;
106 }
107 sbm->blocksize0 = NULL;
108 }
109 if(sbm->blocksize1)
110 {
111 free(sbm->blocksize1);
112 sbm->blocksize1 = NULL;
113 }
114
115 for(unsigned int i = 0 ; i < sbm->nbblocks ; i++)
116 {
117 if(sbm->block[i])
118 {
119 free(sbm->block[i]);
120 sbm->block[i] = NULL;
121 }
122 }
123
124 if(sbm->block)
125 {
126 free(sbm->block);
127 sbm->block = NULL;
128 }
129
130
131 if(sbm->index1_data)
132 {
133 free(sbm->index1_data);
134 sbm->index1_data = NULL;
135 }
136
137 if(sbm->index2_data)
138 {
139 free(sbm->index2_data);
140 sbm->index2_data = NULL;
141 }
142
143 if(sbm->diagonal_blocks)
144 {
145 free(sbm->diagonal_blocks);
146 sbm->diagonal_blocks = NULL;
147 }
148 sbm->filled1 = 0;
149 sbm->filled2 = 0;
150 sbm->blocknumber0 = 0;
151 sbm->blocknumber1 = 0;
152 sbm->nbblocks = 0;
153 }
154
155
SBM_print(const SparseBlockStructuredMatrix * const m)156 void SBM_print(const SparseBlockStructuredMatrix* const m)
157 {
158 if(! m)
159 {
160 fprintf(stderr, "Numerics, SparseBlockStructuredMatrix display failed, NULL input.\n");
161 exit(EXIT_FAILURE);
162 }
163 if(m->blocknumber0 == 0)
164 {
165 printf("Numerics, SparseBlockStructuredMatrix display: matrix dim = 0.");
166 return;
167 }
168
169 int size0 = m->blocksize0[m->blocknumber0 - 1];
170 int size1 = m->blocksize1[m->blocknumber1 - 1];
171 printf("Sparse-Block structured matrix of size %dX%d elements, and %dX%d blocks\n", size0, size1, (int)m->blocknumber0, (int)m->blocknumber1);
172 printf("and %d non null blocks\n", (int)m->nbblocks);
173 printf("Diagonal blocks sizes = [ ");
174 int diagonalblocknumber = m->blocknumber1 + ((m->blocknumber0 - m->blocknumber1) & -(m->blocknumber0 < m->blocknumber1)); // min(m->blocknumber0,m->blocknumber1);
175 for(int i = 0; i < diagonalblocknumber; i++)
176 {
177 size0 = m->blocksize0[i];
178 if(i != 0) size0 -= m->blocksize0[i - 1];
179 size1 = m->blocksize1[i];
180 if(i != 0) size1 -= m->blocksize1[i - 1];
181 printf("%dX%d ", size0, size1);
182 }
183 printf("]\n");
184 printf("index1_data of size %li= {", (long int)m->filled1);
185 if(m->index1_data)
186 {
187 for(unsigned int i = 0 ; i < m->filled1 - 1; i++) printf("%li, ", (long int)m->index1_data[i]);
188 printf("%li}\n", (long int)m->index1_data[m->filled1 - 1]);
189 }
190 else
191 printf("m->index1_data --> NULL}\n");
192
193 printf("index2_data of size %li= {", (long int)m->filled2);
194 if(m->index2_data)
195 {
196 for(unsigned int i = 0 ; i < m->filled2 - 1; i++) printf("%li, ", (long int)m->index2_data[i]);
197 printf("%li}\n", (long int)m->index2_data[m->filled2 - 1]);
198 }
199 else
200 printf("m->index2_data --> NULL}\n");
201
202 printf("blocksize0 of size %li= {", (long int)m->blocknumber0);
203 if(m->blocksize0)
204 {
205 for(unsigned int i = 0 ; i < m->blocknumber0 - 1; i++) printf("%li, ", (long int)m->blocksize0[i]);
206 printf("%li}\n", (long int)m->blocksize0[m->blocknumber0 - 1]);
207 }
208 else
209 printf("m->blocksize0 --> NULL}\n");
210
211 printf("blocksize1 of size %li= {", (long int)m->blocknumber1);
212 if(m->blocksize1)
213 {
214 for(unsigned int i = 0 ; i < m->blocknumber1 - 1; i++) printf("%li, ", (long int)m->blocksize1[i]);
215 printf("%li}\n", (long int)m->blocksize1[m->blocknumber1 - 1]);
216 }
217 else
218 printf("m->blocksize1 --> NULL}\n");
219
220
221 unsigned int currentRowNumber ;
222
223 unsigned int nbRows, nbColumns;
224 if(m->block)
225 {
226 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
227 {
228 /* Get dim. of the current block */
229 nbRows = m->blocksize0[currentRowNumber];
230 if(currentRowNumber != 0)
231 nbRows -= m->blocksize0[currentRowNumber - 1];
232 assert(nbRows);
233
234 for(size_t blockNum = m->index1_data[currentRowNumber];
235 blockNum < m->index1_data[currentRowNumber + 1]; ++blockNum)
236 {
237 assert(blockNum < m->filled2);
238 size_t colNumber = m->index2_data[blockNum];
239 assert(colNumber < m->blocknumber1);
240
241
242
243 nbColumns = m->blocksize1[colNumber];
244 if(colNumber != 0)
245 nbColumns -= m->blocksize1[colNumber - 1];
246 assert(nbColumns);
247
248 printf("block[" SN_SIZE_T_F "] of size %dX%d\n", blockNum, nbRows, nbColumns);
249 if(m->block[blockNum])
250 {
251 unsigned int sizemax = 10;
252 if((nbRows <= sizemax) & (nbColumns <= sizemax))
253 {
254 for(unsigned int i = 0; i < nbRows; i++)
255 {
256 for(unsigned int j = 0; j < nbColumns; j++)
257 {
258 printf("block[" SN_SIZE_T_F "](%i,%i) = %12.8e\n", blockNum, i, j, m->block[blockNum][i + j * nbRows]);
259 }
260 }
261 }
262 else
263 {
264 printf("Block[" SN_SIZE_T_F "] is too large to be displayed\n", blockNum);
265 }
266 }
267 else
268 printf("Block[" SN_SIZE_T_F "] --> NULL \n", blockNum);
269
270 }
271 }
272 }
273 else
274 printf("m->block --> NULL");
275
276 printf("m-diagonal_blocks ");
277 if(m->diagonal_blocks)
278 {
279 printf("[ ");
280 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
281 {
282 printf(" %i",(int)m->diagonal_blocks[currentRowNumber]);
283 }
284 printf("]\n");
285 }
286 else
287 printf("--> NULL");
288
289
290 }
291
292
293 /* a basic iterator scheme for different kind of sparse
294 * matrices (csc, csr, triplet) */
295 typedef struct sparse_matrix_iterator
296 {
297 CS_INT counter1;
298 CS_INT counter2;
299 CS_INT first;
300 CS_INT second;
301 double third;
302 const CSparseMatrix* mat;
303 } sparse_matrix_iterator;
304
305 static sparse_matrix_iterator sparseMatrixBegin(const CSparseMatrix* const sparseMat);
306 static int sparseMatrixNext(sparse_matrix_iterator* it);
307
308
SBM_gemv(unsigned int sizeX,unsigned int sizeY,double alpha,const SparseBlockStructuredMatrix * const restrict A,const double * restrict x,double beta,double * restrict y)309 void SBM_gemv(unsigned int sizeX, unsigned int sizeY, double alpha, const SparseBlockStructuredMatrix* const restrict A, const double* restrict x, double beta, double* restrict y)
310 {
311 /* Product SparseMat - vector, y = A*x (init = 1 = true) or y += A*x (init = 0 = false) */
312
313 assert(A);
314 assert(x);
315 assert(y);
316 assert(A->blocksize0);
317 assert(A->blocksize1);
318 assert(A->index1_data);
319 assert(A->index2_data);
320
321 /* Checks sizes */
322 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
323 assert(sizeY == A->blocksize0[A->blocknumber0 - 1]);
324
325 /* Column (block) position of the current block*/
326 size_t colNumber;
327 /* Number of rows/columns of the current block */
328 unsigned int nbRows, nbColumns;
329 /* Position of the sub-block of x multiplied by the sub-block of A */
330 unsigned int posInX = 0;
331 /* Position of the sub-block of y, result of the product */
332 unsigned int posInY = 0;
333
334 /* Loop over all non-null blocks
335 Works whatever the ordering order of the block is, in A->block
336 */
337 cblas_dscal(sizeY, beta, y, 1);
338
339 for(unsigned int currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
340 {
341 /* Get dim. of the current block */
342 nbRows = A->blocksize0[currentRowNumber];
343 if(currentRowNumber != 0)
344 nbRows -= A->blocksize0[currentRowNumber - 1];
345 assert((nbRows <= sizeY));
346 for(size_t blockNum = A->index1_data[currentRowNumber];
347 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
348 {
349 assert(blockNum < A->filled2);
350
351 colNumber = A->index2_data[blockNum];
352
353 assert(colNumber < sizeX);
354
355 nbColumns = A->blocksize1[colNumber];
356 if(colNumber != 0)
357 nbColumns -= A->blocksize1[colNumber - 1];
358
359 assert((nbColumns <= sizeX));
360
361 /* Get position in x of the sub-block multiplied by A sub-block */
362 posInX = 0;
363 if(colNumber != 0)
364 posInX += A->blocksize1[colNumber - 1];
365 /* Get position in y for the ouput sub-block, result of the product */
366 posInY = 0;
367 if(currentRowNumber != 0)
368 posInY += A->blocksize0[currentRowNumber - 1];
369 /* Computes y[] += currentBlock*x[] */
370 if(nbRows == 3 && nbColumns == 3)
371 {
372 mvp_alpha3x3(alpha, A->block[blockNum], &x[posInX], &y[posInY]);
373 }
374 else
375 {
376 cblas_dgemv(CblasColMajor, CblasNoTrans, nbRows, nbColumns, alpha, A->block[blockNum],
377 nbRows, &x[posInX], 1, 1.0, &y[posInY], 1);
378 }
379 }
380 }
381 }
SBM_gemv_3x3(unsigned int sizeX,unsigned int sizeY,const SparseBlockStructuredMatrix * const restrict A,double * const restrict x,double * restrict y)382 void SBM_gemv_3x3(unsigned int sizeX, unsigned int sizeY, const SparseBlockStructuredMatrix* const restrict A, double* const restrict x, double* restrict y)
383 {
384 /* Product SparseMat - vector, y = vector product y += alpha*A*x for block of size 3x3 */
385
386 assert(A);
387 assert(x);
388 assert(y);
389 assert(A->blocksize0);
390 assert(A->blocksize1);
391 assert(A->index1_data);
392 assert(A->index2_data);
393
394 /* Checks sizes */
395 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
396 assert(sizeY == A->blocksize0[A->blocknumber0 - 1]);
397
398 /* Loop over all non-null blocks
399 Works whatever the ordering order of the block is, in A->block
400 */
401
402 for(unsigned int currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
403 {
404 /* Get dim. of the current block */
405 int nbRows = A->blocksize0[currentRowNumber];
406 if(currentRowNumber != 0)
407 nbRows -= A->blocksize0[currentRowNumber - 1];
408
409 assert((nbRows <= (int)sizeY));
410 for(size_t blockNum = A->index1_data[currentRowNumber];
411 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
412 {
413 assert(blockNum < A->filled2);
414
415 /* Column (block) position of the current block*/
416 size_t colNumber = A->index2_data[blockNum];
417
418 assert(colNumber < sizeX);
419
420 int nbColumns = A->blocksize1[colNumber];
421 if(colNumber != 0)
422 nbColumns -= A->blocksize1[colNumber - 1];
423
424 assert((nbColumns <= (int)sizeX));
425
426 /* Get position in x of the sub-block multiplied by A sub-block */
427 unsigned int posInX = 0;
428 if(colNumber != 0)
429 posInX += A->blocksize1[colNumber - 1];
430
431 /* Get position in y for the ouput sub-block, result of the product */
432 unsigned int posInY = 0;
433 if(currentRowNumber != 0)
434 posInY += A->blocksize0[currentRowNumber - 1];
435
436 /* Computes y[] += currentBlock*x[] */
437
438 /* cblas_dgemv(CblasColMajor, CblasNoTrans, nbRows, nbColumns, alpha, A->block[blockNum], */
439 /* nbRows, &x[posInX], 1, 1.0, &y[posInY], 1); */
440 assert((nbColumns == 3));
441 assert((nbRows == 3));
442 mvp3x3(A->block[blockNum], &x[posInX], &y[posInY]);
443
444 }
445 }
446 }
SBM_extract_component_3x3(const SparseBlockStructuredMatrix * const restrict A,SparseBlockStructuredMatrix * B,unsigned int * row_components,unsigned int row_components_size,unsigned int * col_components,unsigned int col_components_size)447 void SBM_extract_component_3x3(const SparseBlockStructuredMatrix* const restrict A, SparseBlockStructuredMatrix* B,
448 unsigned int *row_components, unsigned int row_components_size,
449 unsigned int *col_components, unsigned int col_components_size)
450 {
451
452
453 assert(A);
454 assert(B);
455 assert(A->blocksize0);
456 assert(A->blocksize1);
457 assert(A->index1_data);
458 assert(A->index2_data);
459
460 /* allocation of data */
461 B->nbblocks= A->nbblocks;
462 B->blocknumber0= A->blocknumber0;
463 B->blocknumber1= A->blocknumber1;
464 B->filled1= A->filled1;
465 B->filled2= A->filled2;
466
467 B->index1_data = (size_t *) malloc(B->filled1*sizeof(size_t));
468 memcpy(B->index1_data, A->index1_data, B->filled1*sizeof(size_t));
469 B->index2_data = (size_t *) malloc(B->filled2*sizeof(size_t));
470 memcpy(B->index2_data, A->index2_data, B->filled2*sizeof(size_t));
471
472 B->blocksize0 = (unsigned int *)malloc(B->blocknumber0*sizeof(unsigned int));
473 B->blocksize1 = (unsigned int *)malloc(B->blocknumber1*sizeof(unsigned int));
474
475
476 int sum =0;
477
478 for(unsigned int row =0; row < B->blocknumber0; row++)
479 {
480 sum +=row_components_size;
481 B->blocksize0[row] = sum ;
482 }
483 sum =0;
484 for(unsigned int col =0; col < B->blocknumber1; col++)
485 {
486 sum +=col_components_size;
487 B->blocksize1[col] = sum ;
488 }
489
490 B->block= (double **)malloc(B->nbblocks*sizeof(double*));
491
492 /* Loop over all non-null blocks
493 Works whatever the ordering order of the block is, in A->block
494 */
495
496 for(unsigned int currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
497 {
498 /* Get dim. of the current block */
499 /* Number of rows of the current block */
500 unsigned int nbRows = A->blocksize0[currentRowNumber];
501 if(currentRowNumber != 0)
502 nbRows -= A->blocksize0[currentRowNumber - 1];
503
504 for(size_t blockNum = A->index1_data[currentRowNumber];
505 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
506 {
507 assert(blockNum < A->filled2);
508
509 B->block[blockNum] = (double*) malloc(row_components_size*col_components_size*sizeof(double));
510
511 for(unsigned int i = 0; i < row_components_size; i++)
512 {
513 for(unsigned int j = 0; j < col_components_size; j++)
514 {
515 B->block[blockNum][i + row_components_size*j] = A->block[blockNum][row_components[i] + col_components[j] * nbRows];
516 }
517 }
518 }
519 }
520
521 DEBUG_EXPR(SBM_print(B));
522 }
523
SBM_check_compatibility_for_add(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B)524 static int SBM_check_compatibility_for_add(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B)
525 {
526 DEBUG_BEGIN(" SBM_check_compatibility_for_add(...)\n");
527 assert(A);
528 assert(B);
529
530 /* Check the compatibility of size of matrices */
531 assert(A->blocksize0);
532 assert(A->blocksize1);
533 assert(B->blocksize0);
534 assert(B->blocksize1);
535
536
537 assert(A->blocksize1[A->blocknumber1 - 1] == B->blocksize1[B->blocknumber1 - 1]);
538 assert(A->blocksize0[A->blocknumber0 - 1] == B->blocksize0[B->blocknumber0 - 1]);
539
540 /* Check the compatibility of the number and the sizes of blocks */
541 int compat = 1;
542
543 if(A->blocknumber0 != B->blocknumber0) compat = 0; /* Compatibility of number of blocks */
544 if(A->blocknumber1 != B->blocknumber1) compat = 0; /* Compatibility of number of blocks */
545
546
547 for(unsigned int i = 0; i < A->blocknumber0; i++)
548 {
549 if(A->blocksize0[i] != B->blocksize0[i]) compat = 0; /* Compatibility of sizes of blocks */
550 }
551 for(unsigned int i = 0; i < A->blocknumber1; i++)
552 {
553 if(A->blocksize1[i] != B->blocksize1[i]) compat = 0; /* Compatibility of sizes of blocks */
554 }
555 DEBUG_PRINTF("compat = %i \n", compat);
556 DEBUG_END("SBM_check_compatibility_for_add(...)\n");
557 return compat;
558
559 }
560
561
562
SBM_calloc_for_add(SparseBlockStructuredMatrix * A,SparseBlockStructuredMatrix * B)563 static SparseBlockStructuredMatrix * SBM_calloc_for_add(SparseBlockStructuredMatrix * A, SparseBlockStructuredMatrix * B)
564 {
565 DEBUG_BEGIN("SBM_calloc_for_add(...)\n");
566 size_t max_number_of_blocks_in_row =0 ;
567
568 for(unsigned int currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
569 {
570
571 size_t nb_blocks_in_rowA = A->index1_data[currentRowNumber + 1] - A->index1_data[currentRowNumber];
572 size_t nb_blocks_in_rowB = B->index1_data[currentRowNumber + 1] - B->index1_data[currentRowNumber];
573 max_number_of_blocks_in_row = max(max_number_of_blocks_in_row, nb_blocks_in_rowA +nb_blocks_in_rowB);
574 }
575
576 size_t * common_block = (size_t *)malloc(max_number_of_blocks_in_row * sizeof(size_t));
577
578 SparseBlockStructuredMatrix * C = SBM_new();
579
580 C->blocknumber0 = A->blocknumber0;
581 C->blocknumber1 = A->blocknumber1;
582 C->blocksize0 = (unsigned int*)malloc(C->blocknumber0*sizeof(unsigned int));
583 C->blocksize1 = (unsigned int*)malloc(C->blocknumber1*sizeof(unsigned int));
584
585 for(unsigned int i = 0; i < A->blocknumber0; i++)
586 {
587 C->blocksize0[i] = A->blocksize0[i];
588 }
589 for(unsigned int i = 0; i < A->blocknumber1; i++)
590 {
591 C->blocksize1[i] = A->blocksize1[i];
592 }
593
594 C->index1_data = (size_t *)malloc(A->filled1 * sizeof(size_t));
595
596 C->index2_data = (size_t *)malloc(max_number_of_blocks_in_row * A->blocknumber0* sizeof(size_t)); /* We oversize a bit */
597 C->block = (double **)malloc(max_number_of_blocks_in_row * A->blocknumber0* sizeof(double *)); /* We oversize a bit */
598
599 C->index1_data[0]=0;
600 C->filled1=0;
601 C->filled2=0;
602
603 for(unsigned int currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
604 {
605
606 size_t nb_blocks_in_rowA = A->index1_data[currentRowNumber + 1] - A->index1_data[currentRowNumber];
607 size_t nb_blocks_in_rowB = B->index1_data[currentRowNumber + 1] - B->index1_data[currentRowNumber];
608 DEBUG_PRINTF("nb_blocks_in_rowA = %zu\n",nb_blocks_in_rowA);
609 DEBUG_PRINTF("nb_blocks_in_rowB = %zu\n",nb_blocks_in_rowB);
610
611 size_t number_of_blocks_in_row = NA_merge_and_sort_sorted_arrays(
612 &(A->index2_data[A->index1_data[currentRowNumber]]),
613 &(B->index2_data[B->index1_data[currentRowNumber]]),nb_blocks_in_rowA,nb_blocks_in_rowB,
614 common_block);
615
616 DEBUG_EXPR(NA_display(common_block,number_of_blocks_in_row););
617
618 DEBUG_PRINTF("number_of_blocks_in_row = %zu\n",number_of_blocks_in_row);
619 C->index1_data[currentRowNumber+1] = C->index1_data[currentRowNumber] + number_of_blocks_in_row ;
620 C->filled1++;
621
622 int blocksize0 = C->blocksize0[currentRowNumber];
623 if(currentRowNumber != 0)
624 blocksize0 -= C->blocksize0[currentRowNumber - 1];
625
626 int idx=0;
627 for(size_t blockNum = C->index1_data[currentRowNumber];
628 blockNum < C->index1_data[currentRowNumber + 1]; ++blockNum)
629 {
630
631 C->filled2++;
632 int currentColNumber = common_block[idx];
633 C->index2_data[blockNum]=currentColNumber;
634 idx++;
635
636 int blocksize1 = C->blocksize1[currentColNumber];
637 if(currentColNumber != 0)
638 blocksize1 -= C->blocksize1[currentColNumber - 1];
639
640 C->block[blockNum] = (double *)calloc(blocksize0*blocksize1,sizeof(double));
641 }
642
643 }
644 C->filled1++;
645 C->nbblocks = C->filled2;
646 DEBUG_END("SBM_calloc_for_add(...)\n");
647 return C;
648 }
649
SBM_add_without_allocation(SparseBlockStructuredMatrix * A,SparseBlockStructuredMatrix * B,double alpha,double beta,SparseBlockStructuredMatrix * C,double gamma)650 void SBM_add_without_allocation(SparseBlockStructuredMatrix * A, SparseBlockStructuredMatrix * B,
651 double alpha, double beta,
652 SparseBlockStructuredMatrix * C, double gamma)
653 {
654 DEBUG_BEGIN("SBM_add_without_allocation(...)\n");
655 assert(SBM_check_compatibility_for_add(A,B) && "Non compatible matrices or blocks sizes.\n");
656
657 size_t colNumber;
658 for(unsigned int currentRowNumber = 0 ; currentRowNumber < C->filled1 - 1; ++currentRowNumber)
659 {
660 DEBUG_PRINT("\n");
661 DEBUG_PRINTF("Computation of the blocks [%i, %i ] in row %i\n", (int)C->index1_data[currentRowNumber], (int)C->index1_data[currentRowNumber+1]-1, (int)currentRowNumber);
662
663 unsigned int blocksize0 = C->blocksize0[currentRowNumber];
664 if(currentRowNumber != 0)
665 blocksize0 -= C->blocksize0[currentRowNumber - 1];
666
667 for(size_t blockNum = C->index1_data[currentRowNumber];
668 blockNum < C->index1_data[currentRowNumber + 1]; ++blockNum)
669 {
670 colNumber=C->index2_data[blockNum];
671 unsigned int blocksize1 = C->blocksize1[colNumber];
672 if(colNumber != 0)
673 blocksize1 -= C->blocksize1[colNumber - 1];
674
675 int nm = blocksize0 * blocksize1;
676 DEBUG_PRINTF("gamma*C for block %zu of size %ix%i\n",blockNum, (int)blocksize0, (int)blocksize1);
677 cblas_dscal(nm, gamma, C->block[blockNum], 1);
678 }
679
680 size_t jC = C->index1_data[currentRowNumber], jA = A->index1_data[currentRowNumber];
681
682 while(jC < C->index1_data[currentRowNumber+1] && jA < A->index1_data[currentRowNumber+1])
683 {
684 /* DEBUG_PRINTF("block number jC = %zu\t, blocknumber jA = %zu\n", jC, jA); */
685 if(C->index2_data[jC] < A->index2_data[jA])
686 {
687 /* DEBUG_PRINT("no block of A to add Block increment jA\n"); */
688 jC++;
689 }
690 else
691 {
692 DEBUG_PRINTF("column number C->index2_data[jC] = %zu\t, column number jA = %zu\n", C->index2_data[jC], A->index2_data[jA]);
693 DEBUG_PRINTF("add a block number %zu of A to blocknumber %zu of C\n",jA,jC);
694
695 unsigned int blocksize1 = C->blocksize1[C->index2_data[jC]];
696 if(C->index2_data[jC] != 0)
697 blocksize1 -= C->blocksize1[C->index2_data[jC] - 1];
698
699 int nm = blocksize0*blocksize1;
700 cblas_daxpy(nm, alpha, A->block[jA], 1, C->block[jC], 1);
701 jA++;
702
703 }
704 }
705 jC = C->index1_data[currentRowNumber];
706 size_t jB = B->index1_data[currentRowNumber];
707
708 while(jC < C->index1_data[currentRowNumber+1] && jB < B->index1_data[currentRowNumber+1])
709 {
710
711 if(C->index2_data[jC] < B->index2_data[jB])
712 {
713 /* DEBUG_PRINT("no block of B to add Block increment jB\n"); */
714 jC++;
715 }
716 else
717 {
718 DEBUG_PRINTF("add a block number %zu of B to blocknumber %zu of C\n",jB,jC);
719 DEBUG_PRINTF("column number C->index2_data[jC] = %zu\t, column number jB = %zu\n", C->index2_data[jC], B->index2_data[jB]);
720 unsigned int blocksize1 = C->blocksize1[C->index2_data[jC]];
721 if(C->index2_data[jC] != 0)
722 blocksize1 -= C->blocksize1[C->index2_data[jC] - 1];
723
724 int nm = blocksize0*blocksize1;
725 cblas_daxpy(nm, beta, B->block[jB], 1, C->block[jC], 1);
726
727 jB++;
728
729 }
730 }
731
732
733
734 }
735 DEBUG_END("SBM_add_without_allocation(...)\n");
736
737 }
SBM_scal(double alpha,SparseBlockStructuredMatrix * A)738 void SBM_scal(double alpha, SparseBlockStructuredMatrix * A)
739 {
740 DEBUG_BEGIN("SBM_scal(...)\n");
741 unsigned int currentRowNumber ;
742 size_t colNumber;
743 unsigned int nbRows, nbColumns;
744
745 for(currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
746 {
747 for(size_t blockNum = A->index1_data[currentRowNumber];
748 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
749 {
750 assert(blockNum < A->filled2);
751 colNumber = A->index2_data[blockNum];
752 /* Get dim. of the current block */
753 nbRows = A->blocksize0[currentRowNumber];
754 if(currentRowNumber != 0)
755 nbRows -= A->blocksize0[currentRowNumber - 1];
756 nbColumns = A->blocksize1[colNumber];
757
758 if(colNumber != 0)
759 nbColumns -= A->blocksize1[colNumber - 1];
760
761 cblas_dscal(nbRows * nbColumns, alpha, A->block[blockNum], 1);
762 }
763 }
764
765 SBM_inc_version(A);
766 DEBUG_END("SBM_scal(...)\n");
767
768 }
769
SBM_add(SparseBlockStructuredMatrix * A,SparseBlockStructuredMatrix * B,double alpha,double beta)770 SparseBlockStructuredMatrix * SBM_add(SparseBlockStructuredMatrix * A, SparseBlockStructuredMatrix * B, double alpha, double beta)
771 {
772
773 /* Check the compatibility of the number and the sizes of blocks */
774 int compat = SBM_check_compatibility_for_add(A,B);
775 if(!compat)
776 {
777 numerics_error("SBM_add", "Non compatible matrices or blocks sizes.\n");
778 return NULL;
779 }
780 SparseBlockStructuredMatrix * C = SBM_calloc_for_add(A,B);
781 SBM_add_without_allocation(A, B, alpha, beta, C, 0.0);
782
783 return C;
784
785 }
786
787
788
789
790
SBM_index_by_column_new()791 static struct SBM_index_by_column* SBM_index_by_column_new()
792 {
793 struct SBM_index_by_column* p = (struct SBM_index_by_column*) malloc(sizeof(struct SBM_index_by_column));
794 p->filled3=0;
795 p->filled4=0;
796 p->index3_data=NULL;
797 p->index4_data=NULL;
798 p->blockMap = NULL;
799 return p;
800 }
SBM_index_by_column_free(struct SBM_index_by_column * p)801 static struct SBM_index_by_column* SBM_index_by_column_free(struct SBM_index_by_column* p)
802 {
803 if(p->index3_data)
804 {
805 free(p->index3_data);
806 p->index3_data=NULL;
807 }
808 if(p->index4_data)
809 {
810 free(p->index4_data);
811 p->index4_data = NULL;
812 }
813 if(p->blockMap)
814 {
815 free(p->blockMap);
816 p->blockMap=NULL;
817 }
818 free(p);
819 p = NULL;
820 return p;
821 }
822
SBM_index_by_column_compute(const SparseBlockStructuredMatrix * const M,struct SBM_index_by_column * SBM_index_by_column_M)823 static void SBM_index_by_column_compute(const SparseBlockStructuredMatrix* const M, struct SBM_index_by_column * SBM_index_by_column_M)
824 {
825
826 SBM_index_by_column_M->filled3 = M->blocknumber1 + 1;
827 SBM_index_by_column_M->index3_data = (size_t *)malloc(SBM_index_by_column_M->filled3 * sizeof(size_t));;
828
829 SBM_index_by_column_M->filled4 = M->nbblocks;
830 SBM_index_by_column_M->index4_data = (size_t *)malloc(SBM_index_by_column_M->filled4 * sizeof(size_t));;
831
832 size_t blockNumM = 0;
833 SBM_index_by_column_M->index3_data[0] = 0;
834 SBM_index_by_column_M->blockMap = (size_t *)malloc(SBM_index_by_column_M->filled4 * sizeof(size_t));
835 unsigned int currentRowNumberofM;
836 unsigned int currentColNumberofM;
837 size_t colNumberofM ;
838 for(currentColNumberofM = 0 ; currentColNumberofM < SBM_index_by_column_M->filled3 - 1; ++currentColNumberofM)
839 {
840 SBM_index_by_column_M->index3_data[currentColNumberofM + 1] = SBM_index_by_column_M->index3_data[currentColNumberofM];
841
842 for(currentRowNumberofM = 0 ; currentRowNumberofM < M->filled1 - 1; ++currentRowNumberofM)
843 {
844 for(size_t blockNum = M->index1_data[currentRowNumberofM];
845 blockNum < M->index1_data[currentRowNumberofM + 1]; ++blockNum)
846 {
847 assert(blockNum < M->filled2);
848 colNumberofM = M->index2_data[blockNum];
849 if(colNumberofM == currentColNumberofM)
850 {
851 assert(blockNumM < M->nbblocks);
852 SBM_index_by_column_M->index3_data[currentColNumberofM + 1]++;
853 SBM_index_by_column_M->index4_data[blockNumM] = currentRowNumberofM;
854 SBM_index_by_column_M->blockMap[blockNumM] = blockNum;
855 blockNumM++;
856 }
857 }
858 }
859 }
860 }
861
SBM_check_compatibility_for_multiply(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B)862 static int SBM_check_compatibility_for_multiply(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B)
863 {
864 DEBUG_BEGIN(" SBM_check_compatibility_for_gemm(...)\n");
865 assert(A);
866 assert(B);
867
868 /* Check the compatibility of size of matrices */
869 assert(A->blocksize0);
870 assert(A->blocksize1);
871 assert(B->blocksize0);
872 assert(B->blocksize1);
873
874 assert(A->blocksize1[A->blocknumber1 - 1] == B->blocksize0[B->blocknumber0 - 1]);
875
876 /* Check the compatibility of the number and the sizes of blocks */
877 int compat = 1;
878
879 if(A->blocknumber1 != B->blocknumber0) compat = 0; /* Compatibility of number of blocks */
880
881
882 for(unsigned int i = 0; i < A->blocknumber1; i++)
883 {
884 if(A->blocksize1[i] != B->blocksize0[i]) compat = 0; /* Compatibility of sizes of blocks */
885 }
886
887 DEBUG_PRINTF("compat = %i \n", compat);
888 DEBUG_END("SBM_check_compatibility_for_gemm(...)\n");
889 return compat;
890
891 }
892
893
SBM_calloc_multiply(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B,struct SBM_index_by_column * SBM_index_by_column_B)894 static SparseBlockStructuredMatrix* SBM_calloc_multiply(const SparseBlockStructuredMatrix* const A,
895 const SparseBlockStructuredMatrix* const B,
896 struct SBM_index_by_column * SBM_index_by_column_B)
897 {
898
899 SparseBlockStructuredMatrix* C = SBM_new();
900
901 C->blocknumber0 = A->blocknumber0;
902 C->blocknumber1 = B->blocknumber1;
903 C->blocksize0 = (unsigned int *)malloc(C->blocknumber0 * sizeof(unsigned int));
904 C->blocksize1 = (unsigned int *)malloc(C->blocknumber1 * sizeof(unsigned int));
905 for(unsigned int i = 0; i < C->blocknumber0; i++) C->blocksize0[i] = A->blocksize0[i];
906 for(unsigned int j = 0; j < C->blocknumber1; j++) C->blocksize1[j] = B->blocksize1[j];
907
908 unsigned int currentRowNumberofA;
909
910 size_t colNumberAA;
911 size_t rowNumberBB;
912 C->nbblocks = 0;
913 C->filled2 = 0;
914 /* \warning The implementation is chosen to optimize cpu effort rather than memory. Otherwise a two loops are needed */
915 int nbblocksmax = A->blocknumber0 * B->blocknumber1;
916
917 double **Cblocktmp = (double**)malloc(nbblocksmax * sizeof(double*));
918 size_t *Cindex2_datatmp = (size_t*)malloc(nbblocksmax * sizeof(size_t));
919 C->filled1 = C->blocknumber0 + 1;
920 C->index1_data = (size_t*)malloc(C->filled1 * sizeof(size_t));
921 C->index1_data[0] = 0;
922 for(currentRowNumberofA = 0 ; currentRowNumberofA < A->filled1 - 1; ++currentRowNumberofA)
923 {
924 C->index1_data[currentRowNumberofA + 1] = C->index1_data[currentRowNumberofA];
925 unsigned int Cblosksize0 = A->blocksize0[currentRowNumberofA];
926 if(currentRowNumberofA != 0)
927 Cblosksize0 -= A->blocksize0[currentRowNumberofA - 1];
928 for(unsigned int currentColNumberofB = 0 ; currentColNumberofB < SBM_index_by_column_B->filled3 - 1; ++currentColNumberofB)
929 {
930 int BlockCexists = 0;
931
932 for(size_t blockNumAA = A->index1_data[currentRowNumberofA];
933 blockNumAA < A->index1_data[currentRowNumberofA + 1]; ++blockNumAA)
934 {
935 assert(blockNumAA < A->filled2);
936 colNumberAA = A->index2_data[blockNumAA];
937
938 for(size_t blockNumBB = SBM_index_by_column_B->index3_data[currentColNumberofB];
939 blockNumBB < SBM_index_by_column_B->index3_data[currentColNumberofB + 1]; ++blockNumBB)
940 {
941 rowNumberBB = SBM_index_by_column_B->index4_data[blockNumBB];
942 if(rowNumberBB == colNumberAA)
943 {
944 BlockCexists = 1;
945 /* printf("C block number %i exists for %i %i ",C->nbblocks, currentRowNumberofA, currentColNumberofB ); */
946
947 unsigned int Cblocksize1 = B->blocksize1[currentColNumberofB];
948 if(currentColNumberofB != 0)
949 Cblocksize1 -= B->blocksize1[currentColNumberofB - 1];
950 /* printf("of size %dX%d\n",Cblosksize0,Cblocksize1 ); */
951 Cblocktmp[C->nbblocks] = (double*)calloc(Cblosksize0 * Cblocksize1, sizeof(double));
952 for(unsigned int i = 0; i < Cblosksize0 * Cblocksize1; i++) Cblocktmp[C->nbblocks][i] = 0.0;
953
954 C->index1_data[currentRowNumberofA + 1]++ ;
955 Cindex2_datatmp[C->nbblocks] = currentColNumberofB ;
956 C->nbblocks++;
957 C->filled2++;
958
959 break;
960
961 };
962 }
963 if(BlockCexists) break;
964 }
965
966 }
967 }
968 assert(C->nbblocks == C->filled2);
969 C->block = (double**)malloc(C->nbblocks * sizeof(double*));
970 C->index2_data = (size_t*)malloc(C->nbblocks * sizeof(size_t));
971
972 for(unsigned int i = 0 ; i < C->nbblocks; i++) C->block[i] = Cblocktmp[i];
973 /* for (i =0 ; i <C->nbblocks; i++) */
974 /* { */
975 /* C->block[i] = (double*)malloc(C->blocksize0[i]*C->blocksize1[i]*sizeof(double)); */
976 /* for (j = 0; j<C->blocksize0[i]*C->blocksize1[i]; j++) C->block[i][j]=0.0; */
977 /* } */
978
979
980 for(unsigned int i = 0 ; i < C->nbblocks; i++) C->index2_data[i] = Cindex2_datatmp[i];
981 free(Cblocktmp);
982 free(Cindex2_datatmp);
983
984 return C;
985 }
SBM_zero_matrix_for_multiply(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B)986 SparseBlockStructuredMatrix* SBM_zero_matrix_for_multiply(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B)
987 {
988 DEBUG_BEGIN("SBM_zero_matrix_for_multiply\n");
989 assert(A);
990 assert(B);
991
992 /* Check the compatibility of the number and the sizes of blocks */
993 int compat = SBM_check_compatibility_for_multiply(A,B);
994 if(!compat)
995 {
996 numerics_error("SBM_zero_matrix_for_multiply", "Non compatible matrices or blocks sizes.\n");
997 return NULL;
998 }
999
1000 /* indexation of B by column */
1001 struct SBM_index_by_column* B_index_by_column = SBM_index_by_column_new();
1002 SBM_index_by_column_compute(B, B_index_by_column);
1003
1004 /* allocation of C */
1005 SparseBlockStructuredMatrix* C = SBM_calloc_multiply(A, B, B_index_by_column);
1006
1007 SBM_index_by_column_free(B_index_by_column);
1008
1009 DEBUG_END("SBM_zero_matrix_for_multiply\n");
1010
1011 return C;
1012 }
1013
1014
1015
SBM_multiply(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B)1016 SparseBlockStructuredMatrix* SBM_multiply(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B)
1017 {
1018 DEBUG_BEGIN("SBM_multiply\(...)\n");
1019 assert(A);
1020 assert(B);
1021
1022
1023 /* Check the compatibility of the number and the sizes of blocks */
1024 int compat = SBM_check_compatibility_for_multiply(A,B);
1025 if(!compat)
1026 {
1027 numerics_error("SBM_multiply", "Non compatible matrices or blocks sizes.\n");
1028 return NULL;
1029 }
1030
1031 /* indexation of B by column */
1032 struct SBM_index_by_column* B_index_by_column = SBM_index_by_column_new();
1033 SBM_index_by_column_compute(B, B_index_by_column);
1034
1035 /* allocation of C */
1036 SparseBlockStructuredMatrix* C = SBM_calloc_multiply(A, B, B_index_by_column);
1037
1038 size_t Bfilled3 = B_index_by_column->filled3;
1039 size_t * Bindex3_data = B_index_by_column->index3_data;
1040 /* unsigned int Bfilled4 = B_index_by_column->filled4; */
1041 size_t * Bindex4_data = B_index_by_column->index4_data;
1042 size_t * blockMap = B_index_by_column->blockMap;
1043
1044 double alpha = 1.0;
1045 double beta = 1.0; /* We assume that C is zero */
1046
1047 size_t colNumberAA;
1048 size_t rowNumberBB;
1049
1050 int Cnbblocks = -1;
1051
1052 for(size_t currentRowNumberofA = 0 ; currentRowNumberofA < A->filled1 - 1; ++currentRowNumberofA)
1053 {
1054 int Ablocksize0 = A->blocksize0[currentRowNumberofA];
1055 if(currentRowNumberofA != 0)
1056 Ablocksize0 -= A->blocksize0[currentRowNumberofA - 1];
1057 for(size_t currentColNumberofB = 0 ; currentColNumberofB < Bfilled3 - 1; ++currentColNumberofB)
1058 {
1059 DEBUG_PRINTF("\n computation of the block of C(%zu,%zu)\n",currentRowNumberofA,currentColNumberofB);
1060 int Bblocksize1 = B->blocksize1[currentColNumberofB];
1061 if(currentColNumberofB != 0)
1062 Bblocksize1 -= B->blocksize1[currentColNumberofB - 1];
1063 /* printf("of size %dX%d\n",Ablocksize0,Bblocksize1 ); */
1064 int CblockPassed = 1;
1065 for(size_t blockNumAA = A->index1_data[currentRowNumberofA];
1066 blockNumAA < A->index1_data[currentRowNumberofA + 1]; ++blockNumAA)
1067 {
1068 assert(blockNumAA < A->filled2);
1069 colNumberAA = A->index2_data[blockNumAA];
1070 DEBUG_PRINTF("blockNumAA = %zu, colNumberAA = %zu\n",blockNumAA,colNumberAA);
1071
1072 for(size_t blockNumBB = Bindex3_data[currentColNumberofB];
1073 blockNumBB < Bindex3_data[currentColNumberofB + 1]; ++blockNumBB)
1074 {
1075 rowNumberBB = Bindex4_data[blockNumBB];
1076 DEBUG_PRINTF("blockNumBB = %zu, rowNumberBB = %zu\n",blockNumBB,rowNumberBB);
1077 DEBUG_PRINTF("blocMap[blockNumBB] = %zu, rowNumberBB = %zu\n",blockMap[blockNumBB],rowNumberBB);
1078
1079 if(rowNumberBB == colNumberAA)
1080 {
1081 if(CblockPassed)
1082 {
1083 Cnbblocks++; /* Find the right C block number*/
1084 CblockPassed = 0;
1085 }
1086 assert(Cnbblocks < (int) C->nbblocks);
1087 /* printf("Compute C block number %i for %i %i ", Cnbblocks,currentRowNumberofA, currentColNumberofB ); */
1088
1089 int Ablocksize1 = A->blocksize1[colNumberAA];
1090 if(colNumberAA != 0)
1091 Ablocksize1 -= A->blocksize1[colNumberAA - 1];
1092
1093 int Bblocksize0 = B->blocksize0[rowNumberBB];
1094 if(rowNumberBB != 0)
1095 Bblocksize0 -= B->blocksize0[rowNumberBB - 1];
1096
1097
1098 DEBUG_PRINTF("Contribution of the product of blocks matrices A(%zu,%zu) (block number %zu) and B(%zu,%zu) (block number %zu) of sizes %dX%d by %dX%d\n",
1099 currentRowNumberofA,colNumberAA, blockNumAA, rowNumberBB, currentColNumberofB, blockMap[blockNumBB],
1100 Ablocksize0,Ablocksize1,Bblocksize0,Bblocksize1);
1101
1102
1103 assert(Ablocksize1 == Bblocksize0);
1104
1105 /* printf("DGEMM call\n"); */
1106 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
1107 Ablocksize0, Bblocksize1, Ablocksize1, alpha, A->block[blockNumAA],
1108 Ablocksize0, B->block[blockMap[blockNumBB]], Bblocksize0,
1109 beta, C->block[Cnbblocks], Ablocksize0);
1110
1111 DEBUG_EXPR(NM_dense_display(A->block[blockNumAA], Ablocksize0, Ablocksize1, 0););
1112 DEBUG_EXPR(NM_dense_display(B->block[blockMap[blockNumBB]], Bblocksize0, Bblocksize1, 0););
1113 DEBUG_EXPR(NM_dense_display(C->block[Cnbblocks], Ablocksize0, Bblocksize1, 0););
1114
1115 /* for (i=0;i<Ablocksize0;i++) */
1116 /* { */
1117 /* for (j=0;j<Bblocksize1;j++) { */
1118 /* printf("C->block[%i](%i,%i) = %f\n",Cnbblocks,i,j, C->block[Cnbblocks][i+j*Ablocksize0]); */
1119 /* } */
1120 /* } */
1121
1122 }
1123 ; /* printf("\n"); */
1124
1125 } /* printf("\n"); */
1126
1127 }
1128
1129 }
1130 }
1131
1132 assert((Cnbblocks + 1) == (int) C->nbblocks);
1133
1134 SBM_index_by_column_free(B_index_by_column);
1135
1136 /* SBM_print(C); */
1137
1138
1139 DEBUG_END("SBM_multiply\(...)\n");
1140 return C;
1141 }
1142
1143 #ifndef NDEBUG
SBM_check_C_for_gemm(const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B,SparseBlockStructuredMatrix * C,struct SBM_index_by_column * SBM_index_by_column_B)1144 static int SBM_check_C_for_gemm(const SparseBlockStructuredMatrix* const A,
1145 const SparseBlockStructuredMatrix* const B,
1146 SparseBlockStructuredMatrix* C,
1147 struct SBM_index_by_column * SBM_index_by_column_B)
1148 {
1149 DEBUG_BEGIN("SBM_check_C_for_gemm(...)\n");
1150
1151 assert(A);
1152 assert(B);
1153 assert(C);
1154
1155 /* Check the compatibility of size of matrices */
1156 assert(A->blocksize0);
1157 assert(A->blocksize1);
1158 assert(B->blocksize0);
1159 assert(B->blocksize1);
1160 assert(C->blocksize0);
1161 assert(C->blocksize1);
1162
1163 assert(A->blocksize1[A->blocknumber1 - 1] == B->blocksize0[B->blocknumber0 - 1]);
1164
1165 /* Check the compatibility of the number and the sizes of blocks */
1166
1167 if(!C)
1168 {
1169 return 0;
1170 }
1171 else
1172 {
1173 if(C->blocknumber0 != A->blocknumber0)
1174 {
1175 DEBUG_PRINT("blocknumber0 are not consistent\n");
1176 return 0;
1177 }
1178 if(C->blocknumber1 != B->blocknumber1)
1179 {
1180 DEBUG_PRINT("blocknumber1 are not consistent\n");
1181 return 0;
1182 }
1183 }
1184
1185 size_t currentRowNumberofA;
1186
1187 size_t colNumberAA;
1188 size_t rowNumberBB;
1189
1190 for(currentRowNumberofA = 0 ; currentRowNumberofA < A->filled1 - 1; ++currentRowNumberofA)
1191 {
1192
1193 if(C->blocksize0[currentRowNumberofA]!= A->blocksize0[currentRowNumberofA])
1194 {
1195 DEBUG_PRINTF("C->blocksize0[%zu] is not consistent with A->blocksize0[%zu]\n",currentRowNumberofA,currentRowNumberofA);
1196 return 0;
1197 }
1198 unsigned int Cblocksize0 = A->blocksize0[currentRowNumberofA];
1199 if(currentRowNumberofA != 0)
1200 Cblocksize0 -= A->blocksize0[currentRowNumberofA - 1];
1201
1202 for(size_t currentColNumberofB = 0 ; currentColNumberofB < SBM_index_by_column_B->filled3 - 1; ++currentColNumberofB)
1203 {
1204 DEBUG_PRINT("\n");
1205 DEBUG_PRINTF("check for the block of C(%zu,%zu)\n",currentRowNumberofA,currentColNumberofB);
1206 int BlockCexists = 0;
1207
1208 if(C->blocksize1[currentColNumberofB]!= B->blocksize1[currentColNumberofB])
1209 {
1210 DEBUG_PRINTF("C->blocksize1[%zu] is not consistent with B->blocksize1[%zu]\n",currentRowNumberofA,currentRowNumberofA);
1211 return 0;
1212 }
1213 for(size_t blockNumAA = A->index1_data[currentRowNumberofA];
1214 blockNumAA < A->index1_data[currentRowNumberofA + 1]; ++blockNumAA)
1215 {
1216 assert(blockNumAA < A->filled2);
1217 colNumberAA = A->index2_data[blockNumAA];
1218
1219 for(size_t blockNumBB = SBM_index_by_column_B->index3_data[currentColNumberofB];
1220 blockNumBB < SBM_index_by_column_B->index3_data[currentColNumberofB + 1]; ++blockNumBB)
1221 {
1222 rowNumberBB = SBM_index_by_column_B->index4_data[blockNumBB];
1223
1224 if(rowNumberBB == colNumberAA)
1225 {
1226 BlockCexists = 1;
1227
1228 unsigned int Cblocksize1 = B->blocksize1[currentColNumberofB];
1229 if(currentColNumberofB != 0)
1230 Cblocksize1 -= B->blocksize1[currentColNumberofB - 1];
1231 DEBUG_PRINTF("C block number is needed for %zu %zu of size %dX%d\n",
1232 currentRowNumberofA, currentColNumberofB, (int)Cblocksize0, (int)Cblocksize1);
1233
1234 /* search for the block */
1235 size_t C_blockNum = 0;
1236 int block_C_found = 0;
1237 /* find the correct blockNum is C */
1238 for(size_t blockNumCC= C->index1_data[currentRowNumberofA];
1239 blockNumCC < C->index1_data[currentRowNumberofA + 1]; ++blockNumCC)
1240 {
1241 if(C->index2_data[blockNumCC] == currentColNumberofB)
1242 {
1243 C_blockNum = blockNumCC;
1244 block_C_found = 1;
1245 DEBUG_PRINTF("Block in C (number %zu ) is found with corresponding row %zu and col %zu\n",
1246 C_blockNum, currentRowNumberofA, currentColNumberofB);
1247 }
1248 }
1249 if(!block_C_found)
1250 {
1251 DEBUG_PRINTF("Block in C (number %zu ) is found with corresponding row %zu and col %zu\n",
1252 C_blockNum, currentRowNumberofA, currentColNumberofB);
1253 return 0;
1254 }
1255
1256 if(!C->block[C_blockNum])
1257 {
1258 DEBUG_PRINTF("Block in C (number %zu ) is found but not allocated\n",
1259 C_blockNum);
1260 return 0;
1261 }
1262
1263 break;
1264 };
1265 }
1266 if(BlockCexists) break;
1267 }
1268 }
1269 }
1270
1271 DEBUG_END("SBM_check_C_for_gemm(...)\n");
1272 return 1;
1273 }
1274 #endif
1275
SBM_gemm_without_allocation(double alpha,const SparseBlockStructuredMatrix * const A,const SparseBlockStructuredMatrix * const B,double beta,SparseBlockStructuredMatrix * C)1276 void SBM_gemm_without_allocation(double alpha, const SparseBlockStructuredMatrix* const A,
1277 const SparseBlockStructuredMatrix* const B,
1278 double beta, SparseBlockStructuredMatrix* C)
1279 {
1280 DEBUG_BEGIN("SBM_gemm_without_allocation\(...)\n");
1281
1282 /* Check the compatibility of the number and the sizes of blocks */
1283 assert(SBM_check_compatibility_for_multiply(A,B));
1284
1285 struct SBM_index_by_column* B_index_by_column = SBM_index_by_column_new();
1286 SBM_index_by_column_compute(B, B_index_by_column);
1287
1288
1289 /* Check the memory allocation in C */
1290 assert(SBM_check_C_for_gemm(A, B, C, B_index_by_column));
1291
1292 size_t Bfilled3 = B_index_by_column->filled3;
1293 size_t * Bindex3_data = B_index_by_column->index3_data;
1294 /* unsigned int Bfilled4 = B_index_by_column->filled4; */
1295 size_t * Bindex4_data = B_index_by_column->index4_data;
1296 size_t * blockMap = B_index_by_column->blockMap;
1297
1298 size_t colNumberAA;
1299 size_t rowNumberBB;
1300
1301 size_t C_blockNum = SIZE_MAX;
1302 for(size_t currentRowNumberofA = 0 ; currentRowNumberofA < A->filled1 - 1; ++currentRowNumberofA)
1303 {
1304 int Ablocksize0 = A->blocksize0[currentRowNumberofA];
1305 if(currentRowNumberofA != 0)
1306 Ablocksize0 -= A->blocksize0[currentRowNumberofA - 1];
1307 for(size_t currentColNumberofB = 0 ; currentColNumberofB < Bfilled3 - 1; ++currentColNumberofB)
1308 {
1309 DEBUG_PRINTF("\n computation of the block of C(%zu,%zu)\n",currentRowNumberofA,currentColNumberofB);
1310
1311 int Bblocksize1 = B->blocksize1[currentColNumberofB];
1312 if(currentColNumberofB != 0)
1313 Bblocksize1 -= B->blocksize1[currentColNumberofB - 1];
1314 int CblockPassed = 1;
1315 double local_beta=beta;
1316 for(size_t blockNumAA = A->index1_data[currentRowNumberofA];
1317 blockNumAA < A->index1_data[currentRowNumberofA + 1]; ++blockNumAA)
1318 {
1319 assert(blockNumAA < A->filled2);
1320 colNumberAA = A->index2_data[blockNumAA];
1321 DEBUG_PRINTF("blockNumAA = %zu, colNumberAA = %zu\n",blockNumAA,colNumberAA);
1322
1323
1324 for(size_t blockNumBB = Bindex3_data[currentColNumberofB];
1325 blockNumBB < Bindex3_data[currentColNumberofB + 1]; ++blockNumBB)
1326 {
1327 rowNumberBB = Bindex4_data[blockNumBB];
1328 DEBUG_PRINTF("blockNumBB = %zu, rowNumberBB = %zu\n",blockNumBB,rowNumberBB);
1329 DEBUG_PRINTF("blocMap[blockNumBB] = %zu, rowNumberBB = %zu\n",blockMap[blockNumBB],rowNumberBB);
1330
1331
1332 if(rowNumberBB == colNumberAA)
1333 {
1334 /* search for the block */
1335 if(CblockPassed)
1336 {
1337 CblockPassed = 0;
1338 C_blockNum = SIZE_MAX;
1339 /* find the correct blockNum is C */
1340 for(size_t blockNumCC= C->index1_data[currentRowNumberofA];
1341 blockNumCC < C->index1_data[currentRowNumberofA + 1]; ++blockNumCC)
1342 {
1343 if(C->index2_data[blockNumCC] == currentColNumberofB)
1344 {
1345 C_blockNum = blockNumCC;
1346 DEBUG_PRINTF("Block in C (number %zu ) is found with corresponding row %zu and col %zu\n",
1347 C_blockNum, currentRowNumberofA, currentColNumberofB);
1348 }
1349 }
1350 }
1351 assert(C_blockNum != SIZE_MAX);
1352
1353
1354
1355 /* printf("of size %dX%d\n",Ablocksize0,Bblocksize1 ); */
1356
1357 int Ablocksize1 = A->blocksize1[colNumberAA];
1358 if(colNumberAA != 0)
1359 Ablocksize1 -= A->blocksize1[colNumberAA - 1];
1360
1361 int Bblocksize0 = B->blocksize0[rowNumberBB];
1362 if(rowNumberBB != 0)
1363 Bblocksize0 -= B->blocksize0[rowNumberBB - 1];
1364
1365
1366 DEBUG_PRINTF("Contribution of the product of blocks matrices A(%zu,%zu) (block number %zu) and B(%zu,%zu) (block number %zu) of sizes %dX%d by %dX%d\n",
1367 currentRowNumberofA,colNumberAA, blockNumAA, rowNumberBB, currentColNumberofB, blockMap[blockNumBB],
1368 Ablocksize0,Ablocksize1,Bblocksize0,Bblocksize1);
1369
1370 assert(Ablocksize1 == Bblocksize0);
1371
1372 /* printf("DGEMM call\n"); */
1373 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
1374 Ablocksize0, Bblocksize1, Ablocksize1, alpha, A->block[blockNumAA],
1375 Ablocksize0, B->block[blockMap[blockNumBB]], Bblocksize0,
1376 local_beta, C->block[C_blockNum], Ablocksize0);
1377 local_beta=1.0;
1378 DEBUG_EXPR(NM_dense_display(A->block[blockNumAA], Ablocksize0, Ablocksize1, 0););
1379 DEBUG_EXPR(NM_dense_display(B->block[blockMap[blockNumBB]], Bblocksize0, Bblocksize1, 0););
1380 DEBUG_EXPR(NM_dense_display(C->block[C_blockNum], Ablocksize0, Bblocksize1, 0););
1381 }
1382 }
1383 }
1384 }
1385 }
1386
1387 SBM_index_by_column_free(B_index_by_column);
1388
1389 DEBUG_END("SBM_gemm_without_allocation\(...)\n")
1390 }
SBM_row_prod(unsigned int sizeX,unsigned int sizeY,unsigned int currentRowNumber,const SparseBlockStructuredMatrix * const A,const double * const x,double * y,int init)1391 void SBM_row_prod(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber,
1392 const SparseBlockStructuredMatrix* const A,
1393 const double* const x, double* y, int init)
1394 {
1395 /*
1396 Product (Row of blocks of a SparseMat) - vector, y = rowA*x (init
1397 = 1 = true) or y += rowA*x (init = 0 = false)
1398 */
1399
1400
1401 assert(A);
1402 assert(x);
1403 assert(y);
1404
1405 /* Checks sizes */
1406 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
1407
1408 /* Column (block) position of the current block*/
1409 size_t colNumber = 0;
1410
1411 /* Check if currentRowNumber fits with A dimensions */
1412 assert(currentRowNumber <= A->blocknumber0);
1413
1414 /* Get dim (rows) of the current block */
1415 unsigned int nbRows = sizeY;
1416
1417 /* if this is important, move it into a function --xhub */
1418 /*
1419 assert(
1420 {
1421 nbRows = A->blocksize0[currentRowNumber];
1422 if(currentRowNumber!=0)
1423 nbRows -= A->blocksize0[currentRowNumber-1];
1424 nbRows ==sizeY;
1425 }); */
1426
1427 /* Set y to 0, if required */
1428 if(init == 1)
1429 cblas_dscal(sizeY, 0.0, y, 1);
1430
1431 /* Loop over all non-null blocks
1432 Works whatever the ordering order of the block is, in A->block
1433 But it requires a set to 0 of all y components
1434 */
1435
1436 for(size_t blockNum = A->index1_data[currentRowNumber];
1437 blockNum < A->index1_data[currentRowNumber + 1];
1438 ++blockNum)
1439 {
1440 /* Get row/column position of the current block */
1441 colNumber = A->index2_data[blockNum];
1442
1443 /* Get dim(columns) of the current block */
1444 unsigned int nbColumns = A->blocksize1[colNumber];
1445 if(colNumber != 0)
1446 nbColumns -= A->blocksize1[colNumber - 1];
1447
1448 /* Get position in x of the sub-block multiplied by A sub-block */
1449 /* Position of the sub-block of x multiplied by the sub-block of A */
1450 unsigned int posInX =0;
1451 if(colNumber != 0)
1452 posInX += A->blocksize0[colNumber - 1];
1453 /* Computes y[] += currentBlock*x[] */
1454 cblas_dgemv(CblasColMajor, CblasNoTrans, nbRows, nbColumns, 1.0, A->block[blockNum], nbRows, &x[posInX], 1, 1.0, y, 1);
1455
1456 }
1457 }
SBM_row_prod_no_diag(unsigned int sizeX,unsigned int sizeY,unsigned int currentRowNumber,const SparseBlockStructuredMatrix * const A,const double * const x,double * y,int init)1458 void SBM_row_prod_no_diag(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, const double* const x, double* y, int init)
1459 {
1460 /*
1461 If: A is a SparseBlockStructuredMatrix matrix, Aij a block at row
1462 i and column j (Warning: i and j are indices of block position,
1463 not scalar component positions)
1464
1465 Then SBM_row_prod_no_diag computes y = sum for i not equal to j of
1466 Aij.xj over a row of blocks (or += if init = false)
1467
1468 currentRowNumber represents the position (block number) of the
1469 required line of blocks in the matrix A.
1470
1471 */
1472
1473
1474 /* Column (block) position of the current block*/
1475 size_t colNumber = 0;
1476
1477 /* Number of rows/columns of the current block */
1478 unsigned int nbRows, nbColumns;
1479
1480 /* Position of the sub-block of x multiplied by the sub-block of
1481 * A */
1482 unsigned int posInX = 0;
1483
1484 /* Look for the first element of the wanted row */
1485
1486 /* Assertions */
1487 assert(A);
1488 assert(x);
1489 assert(y);
1490 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
1491 assert(currentRowNumber <= A->blocknumber0);
1492
1493 /* Get dim (rows) of the current block */
1494 nbRows = sizeY;
1495
1496 /* if this is important, move it into a function --xhub */
1497 /* assert(
1498 {
1499 nbRows = A->blocksize0[currentRowNumber];
1500 if(currentRowNumber!=0)
1501 nbRows -= A->blocksize0[currentRowNumber-1];
1502 nbRows == sizeY ;
1503 });*/
1504
1505 /* Set y to 0, if required */
1506 if(init == 1)
1507 cblas_dscal(sizeY, 0.0, y, 1);
1508
1509 /* Loop over all non-null blocks. Works whatever the ordering order
1510 of the block is, in A->block, but it requires a set to 0 of all y
1511 components
1512 */
1513 for(size_t blockNum = A->index1_data[currentRowNumber];
1514 blockNum < A->index1_data[currentRowNumber + 1];
1515 ++blockNum)
1516 {
1517 /* Get row/column position of the current block */
1518 colNumber = A->index2_data[blockNum];
1519
1520 /* Computes product only for extra diagonal blocks */
1521 if(colNumber != currentRowNumber)
1522 {
1523 /* Get dim(columns) of the current block */
1524 nbColumns = A->blocksize1[colNumber];
1525 if(colNumber != 0)
1526 nbColumns -= A->blocksize1[colNumber - 1];
1527
1528 /* Get position in x of the sub-block multiplied by A sub-block */
1529 posInX = 0;
1530 if(colNumber != 0)
1531 posInX += A->blocksize0[colNumber - 1];
1532 /* Computes y[] += currentBlock*x[] */
1533 cblas_dgemv(CblasColMajor,CblasNoTrans, nbRows, nbColumns, 1.0, A->block[blockNum], nbRows, &x[posInX], 1, 1.0, y, 1);
1534 }
1535 }
1536 }
SBM_row_prod_no_diag_3x3(unsigned int sizeX,unsigned int sizeY,unsigned int currentRowNumber,const SparseBlockStructuredMatrix * const A,double * const x,double * y)1537 void SBM_row_prod_no_diag_3x3(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, double* const x, double* y)
1538 {
1539 /*
1540 If: A is a SparseBlockStructuredMatrix matrix, Aij a block at row
1541 i and column j (Warning: i and j are indices of block position,
1542 not scalar component positions)
1543
1544 Then SBM_row_prod_no_diag computes y = sum for i not equal to j of
1545 Aij.xj over a row of blocks (or += if init = false)
1546
1547 currentRowNumber represents the position (block number) of the
1548 required line of blocks in the matrix A.
1549 */
1550
1551 /* Number of columns of the current block */
1552 unsigned int nbColumns;
1553
1554 /* Position of the sub-block of x multiplied by the sub-block of
1555 * A */
1556 unsigned int posInX = 0;
1557
1558 /* Look for the first element of the wanted row */
1559
1560 /* Assertions */
1561 assert(A);
1562 assert(x);
1563 assert(y);
1564 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
1565 assert(currentRowNumber <= A->blocknumber0);
1566
1567 /* Loop over all non-null blocks. Works whatever the ordering order
1568 of the block is, in A->block, but it requires a set to 0 of all y
1569 components
1570 */
1571 for(size_t blockNum = A->index1_data[currentRowNumber];
1572 blockNum < A->index1_data[currentRowNumber + 1];
1573 ++blockNum)
1574 {
1575 /* Get row/column position of the current block */
1576 size_t colNumber = A->index2_data[blockNum];
1577
1578 /* Computes product only for extra diagonal blocks */
1579 if(colNumber != currentRowNumber)
1580 {
1581 /* Get dim(columns) of the current block */
1582 nbColumns = A->blocksize1[colNumber];
1583 if(colNumber != 0)
1584 nbColumns -= A->blocksize1[colNumber - 1];
1585
1586 /* Get position in x of the sub-block multiplied by A sub-block */
1587 posInX = 0;
1588 if(colNumber != 0)
1589 posInX += A->blocksize0[colNumber - 1];
1590 /* Computes y[] += currentBlock*x[] */
1591 /* cblas_dgemv(CblasColMajor,CblasNoTrans, nbRows, nbColumns, 1.0, A->block[blockNum], nbRows, &x[posInX], 1, 1.0, y, 1); */
1592 assert((nbColumns == 3));
1593 mvp3x3(A->block[blockNum], &x[posInX], y);
1594 }
1595 }
1596 }
SBM_row_prod_no_diag_1x1(unsigned int sizeX,unsigned int sizeY,unsigned int currentRowNumber,const SparseBlockStructuredMatrix * const A,double * const x,double * y)1597 void SBM_row_prod_no_diag_1x1(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, double* const x, double* y)
1598 {
1599 /*
1600 If: A is a SparseBlockStructuredMatrix matrix, Aij a block at row
1601 i and column j (Warning: i and j are indices of block position,
1602 not scalar component positions)
1603
1604 Then SBM_row_prod_no_diag computes y = sum for i not equal to j of
1605 Aij.xj over a row of blocks (or += if init = false)
1606
1607 currentRowNumber represents the position (block number) of the
1608 required line of blocks in the matrix A.
1609
1610 */
1611
1612 /* Number of columns of the current block */
1613 unsigned int nbColumns;
1614
1615
1616 /* Look for the first element of the wanted row */
1617
1618 /* Assertions */
1619 assert(A);
1620 assert(x);
1621 assert(y);
1622 assert(sizeX == A->blocksize1[A->blocknumber1 - 1]);
1623 assert(currentRowNumber <= A->blocknumber0);
1624
1625 /* Loop over all non-null blocks. Works whatever the ordering order
1626 of the block is, in A->block, but it requires a set to 0 of all y
1627 components
1628 */
1629 for(size_t blockNum = A->index1_data[currentRowNumber];
1630 blockNum < A->index1_data[currentRowNumber + 1];
1631 ++blockNum)
1632 {
1633 /* Get row/column position of the current block */
1634 size_t colNumber = A->index2_data[blockNum];
1635
1636 /* Computes product only for extra diagonal blocks */
1637 if(colNumber != currentRowNumber)
1638 {
1639 /* Get dim(columns) of the current block */
1640 nbColumns = A->blocksize1[colNumber];
1641 if(colNumber != 0)
1642 nbColumns -= A->blocksize1[colNumber - 1];
1643
1644 /* Get position in x of the sub-block multiplied by A sub-block */
1645 unsigned int posInX = 0;
1646 if(colNumber != 0)
1647 posInX += A->blocksize0[colNumber - 1];
1648 /* Computes y[] += currentBlock*x[] */
1649 /* cblas_dgemv(CblasColMajor,CblasNoTrans, nbRows, nbColumns, 1.0, A->block[blockNum], nbRows, &x[posInX], 1, 1.0, y, 1); */
1650 assert((nbColumns == 1));
1651 y[0] += A->block[blockNum][0] * x[posInX] ;
1652 }
1653 }
1654 }
1655
1656
SBM_write_in_file(const SparseBlockStructuredMatrix * const m,FILE * file)1657 void SBM_write_in_file(const SparseBlockStructuredMatrix* const m, FILE * file)
1658 {
1659 DEBUG_PRINT("printInFileSBM\n");
1660 if(! m)
1661 {
1662 fprintf(stderr, "Numerics, SparseBlockStructuredMatrix printInFileSBM failed, NULL input.\n");
1663 exit(EXIT_FAILURE);
1664 }
1665 assert(file);
1666 fprintf(file, "%i\n", (int)m->nbblocks);
1667 if(m->nbblocks == 0) return;
1668
1669 fprintf(file, "%i\n", (int)m->blocknumber0);
1670 fprintf(file, "%i\n", (int)m->blocknumber1);
1671
1672 assert(m->blocksize0);
1673 assert(m->blocksize1);
1674 assert(m->index1_data);
1675 assert(m->index2_data);
1676
1677 for(unsigned int i = 0; i < m->blocknumber0; i++)
1678 {
1679 fprintf(file, "%i\t", m->blocksize0[i]);
1680 }
1681 fprintf(file, "\n");
1682 for(unsigned int i = 0; i < m->blocknumber1; i++)
1683 {
1684 fprintf(file, "%i\t", m->blocksize1[i]);
1685 }
1686 fprintf(file, "\n");
1687 fprintf(file, "%li\n", (long int)m->filled1);
1688 fprintf(file, "%li\n", (long int)m->filled2);
1689 for(unsigned int i = 0; i < m->filled1; i++)
1690 {
1691 fprintf(file, "%li\t", (long int)m->index1_data[i]);
1692 }
1693 fprintf(file, "\n");
1694 for(unsigned int i = 0; i < m->filled2; i++)
1695 {
1696 fprintf(file, "%li\t", (long int)m->index2_data[i]);
1697 }
1698 fprintf(file, "\n");
1699
1700 unsigned int currentRowNumber ;
1701 size_t colNumber;
1702 unsigned int nbRows, nbColumns;
1703 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
1704 {
1705 DEBUG_PRINTF("currentRowNumber = %i\n", (int)currentRowNumber);
1706 DEBUG_PRINTF(" m->index1_data[currentRowNumber] = %zu\n", m->index1_data[currentRowNumber]);
1707 DEBUG_PRINTF(" m->index1_data[currentRowNumber+1] = %zu\n", m->index1_data[currentRowNumber+1]);
1708
1709 for(size_t blockNum = m->index1_data[currentRowNumber];
1710 blockNum < m->index1_data[currentRowNumber + 1]; ++blockNum)
1711 {
1712 assert(blockNum < m->filled2);
1713 colNumber = m->index2_data[blockNum];
1714 /* Get dim. of the current block */
1715 nbRows = m->blocksize0[currentRowNumber];
1716
1717 if(currentRowNumber != 0)
1718 nbRows -= m->blocksize0[currentRowNumber - 1];
1719
1720 nbColumns = m->blocksize1[colNumber];
1721 if(colNumber != 0)
1722 nbColumns -= m->blocksize1[colNumber - 1];
1723 //fprintf(file,"block[%i] of size %dX%d\n", blockNum, nbRows,nbColumns);
1724 fprintf(file, SN_SIZE_T_F "\n", blockNum);
1725 DEBUG_PRINTF("nbRows * nbColumns = %i\n", (int)(nbRows * nbColumns));
1726 assert(m->block[blockNum]);
1727 for(unsigned int i = 0; i < nbRows * nbColumns; i++)
1728 {
1729 DEBUG_PRINTF("i = %i, blockNum = %zu, m->block[%zu][%i] = %g\n ", (int)i, blockNum, blockNum, i, m->block[blockNum][i]);
1730 fprintf(file, "%32.24e\n", m->block[blockNum][i]);
1731 }
1732
1733 }
1734 }
1735
1736 }
SBM_write_in_fileForScilab(const SparseBlockStructuredMatrix * const m,FILE * file)1737 void SBM_write_in_fileForScilab(const SparseBlockStructuredMatrix* const m, FILE * file)
1738 {
1739 if(! m)
1740 {
1741 fprintf(stderr, "Numerics, SparseBlockStructuredMatrix SBM_write_in_file failed, NULL input.\n");
1742 exit(EXIT_FAILURE);
1743 }
1744 assert(file);
1745 fprintf(file, "nbblock = %i;\n", (int)m->nbblocks);
1746 if(m->nbblocks == 0) return;
1747
1748 fprintf(file, "blocknumber0 = %i;\n", (int)m->blocknumber0);
1749 fprintf(file, "blocknumber1 = %i; \n", (int)m->blocknumber1);
1750
1751 assert(m->blocksize0);
1752 assert(m->blocksize1);
1753 assert(m->index1_data);
1754 assert(m->index2_data);
1755
1756 fprintf(file, "blocksize0 = [ \t");
1757 for(unsigned int i = 0; i < m->blocknumber0; i++)
1758 {
1759 fprintf(file, "%i\t", (int)m->blocksize0[i]);
1760 }
1761 fprintf(file, "];\n");
1762 fprintf(file, "blocksize1 = [ \t");
1763 for(unsigned int i = 0; i < m->blocknumber1; i++)
1764 {
1765 fprintf(file, "%i\t", (int)m->blocksize1[i]);
1766 }
1767 fprintf(file, "];\n");
1768 fprintf(file, "filled1 = %li;\n", (long int)m->filled1);
1769 fprintf(file, "filled2 = %li;\n", (long int)m->filled2);
1770
1771 fprintf(file, "index1_data = [ \t");
1772 for(unsigned int i = 0; i < m->filled1; i++)
1773 {
1774 fprintf(file, "%li\t", (long int)m->index1_data[i]);
1775 }
1776 fprintf(file, "];\n");
1777 fprintf(file, "index2_data = [ \t");
1778 for(unsigned int i = 0; i < m->filled2; i++)
1779 {
1780 fprintf(file, "%li\t", (long int)m->index2_data[i]);
1781 }
1782 fprintf(file, "];\n");
1783
1784 unsigned int currentRowNumber ;
1785 size_t colNumber;
1786 unsigned int nbRows, nbColumns;
1787 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
1788 {
1789 for(size_t blockNum = m->index1_data[currentRowNumber];
1790 blockNum < m->index1_data[currentRowNumber + 1]; ++blockNum)
1791 {
1792 assert(blockNum < m->filled2);
1793 colNumber = m->index2_data[blockNum];
1794 /* Get dim. of the current block */
1795 nbRows = m->blocksize0[currentRowNumber];
1796
1797 if(currentRowNumber != 0)
1798 nbRows -= m->blocksize0[currentRowNumber - 1];
1799
1800 nbColumns = m->blocksize1[colNumber];
1801 if(colNumber != 0)
1802 nbColumns -= m->blocksize1[colNumber - 1];
1803 //fprintf(file,"block[%i] of size %dX%d\n", blockNum, nbRows,nbColumns);
1804 fprintf(file, "block" SN_SIZE_T_F " = [ \n", blockNum);
1805
1806 for(unsigned int i = 0; i < nbRows; i++)
1807 {
1808 fprintf(file, "[");
1809 for(unsigned int j = 0; j < nbColumns; j++)
1810 {
1811 fprintf(file, "%32.24e\t ", m->block[blockNum][i + j * nbRows]);
1812 }
1813 fprintf(file, "];\n");
1814 }
1815 fprintf(file, "];\n");
1816 /* for (int i=0;i<nbRows*nbColumns;i++) */
1817 /* { */
1818 /* fprintf(file,"%32le\n",m->block[blockNum][i]); */
1819 /* } */
1820 }
1821 }
1822
1823 fprintf(file, "// Dense version\n");
1824 int size0 = m->blocksize0[m->blocknumber0 - 1];
1825 int size1 = m->blocksize1[m->blocknumber1 - 1];
1826 double * denseMat = (double*)malloc(size0 * size1 * sizeof(double));
1827 SBM_to_dense(m, denseMat);
1828
1829 fprintf(file, "data= [");
1830 for(int i = 0; i < size0; i++)
1831 {
1832 fprintf(file, "[");
1833 for(int j = 0; j < size1; j++)
1834 {
1835 fprintf(file, "%32.24e,\t ", denseMat[i + j * size0]);
1836 }
1837 fprintf(file, "];\n");
1838 }
1839 fprintf(file, "]");
1840 free(denseMat);
1841
1842
1843
1844 }
1845
1846
SBM_write_in_filename(const SparseBlockStructuredMatrix * const m,const char * filename)1847 void SBM_write_in_filename(const SparseBlockStructuredMatrix* const m, const char *filename)
1848 {
1849
1850 }
SBM_new_from_file(FILE * file)1851 SparseBlockStructuredMatrix * SBM_new_from_file(FILE *file)
1852 {
1853
1854 SparseBlockStructuredMatrix * m = SBM_new();
1855
1856 assert(file);
1857 CHECK_IO(fscanf(file, "%d", &(m->nbblocks)));
1858
1859 if(m->nbblocks == 0) return NULL;
1860
1861 CHECK_IO(fscanf(file, "%d", &(m->blocknumber0)));
1862 CHECK_IO(fscanf(file, "%d", &(m->blocknumber1)));
1863
1864 m->blocksize0 = (unsigned int *)malloc(m->blocknumber0 * sizeof(unsigned int));
1865 m->blocksize1 = (unsigned int *)malloc(m->blocknumber1 * sizeof(unsigned int));
1866
1867
1868 for(unsigned int i = 0; i < m->blocknumber0; i++)
1869 {
1870 CHECK_IO(fscanf(file, "%d", &(m->blocksize0[i])));
1871 }
1872 for(unsigned int i = 0; i < m->blocknumber1; i++)
1873 {
1874 CHECK_IO(fscanf(file, "%d", &(m->blocksize1[i])));
1875 }
1876
1877 unsigned int filled1 = 0, filled2 = 0;
1878 CHECK_IO(fscanf(file, "%d", &(filled1)));
1879 m->filled1 = filled1;
1880 CHECK_IO(fscanf(file, "%d", &(filled2)));
1881 m->filled2 = filled2;
1882 m->index1_data = (size_t *)malloc(m->filled1 * sizeof(size_t));
1883 m->index2_data = (size_t *)malloc(m->filled2 * sizeof(size_t));
1884
1885 int index1_dataCurrent = 0;
1886 for(unsigned int i = 0; i < m->filled1; i++)
1887 {
1888 CHECK_IO(fscanf(file, "%d", &(index1_dataCurrent)));
1889 m->index1_data[i] = index1_dataCurrent;
1890 }
1891 int index2_dataCurrent = 0;
1892 for(unsigned int i = 0; i < m->filled2; i++)
1893 {
1894 CHECK_IO(fscanf(file, "%d", &(index2_dataCurrent)));
1895 m->index2_data[i] = index2_dataCurrent;
1896 }
1897 m->block = (double**)malloc(m->nbblocks * sizeof(double*));
1898 unsigned int currentRowNumber ;
1899 size_t colNumber;
1900 unsigned int nbRows, nbColumns;
1901 unsigned int blockk;
1902 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
1903 {
1904 for(size_t blockNum = m->index1_data[currentRowNumber];
1905 blockNum < m->index1_data[currentRowNumber + 1]; ++blockNum)
1906 {
1907 assert(blockNum < m->filled2);
1908 colNumber = m->index2_data[blockNum];
1909 /* Get dim. of the current block */
1910 nbRows = m->blocksize0[currentRowNumber];
1911
1912 if(currentRowNumber != 0)
1913 nbRows -= m->blocksize0[currentRowNumber - 1];
1914
1915 nbColumns = m->blocksize1[colNumber];
1916 if(colNumber != 0)
1917 nbColumns -= m->blocksize1[colNumber - 1];
1918 //fprintf(file,"block[%i] of size %dX%d\n", blockNum, nbRows,nbColumns);
1919
1920 CHECK_IO(fscanf(file, "%d", &(blockk)));
1921 if(blockk != blockNum)
1922 {
1923 printf("Numerics, SparseBlockStructuredMatrix SBM_read_in_file failed, problem in block numbering. \n");
1924 }
1925 m->block[blockNum] = (double*)malloc(nbRows * nbColumns * sizeof(double));
1926 for(unsigned int i = 0; i < nbRows * nbColumns; i++)
1927 {
1928 CHECK_IO(fscanf(file, "%32le\n", &(m->block[blockNum][i])));
1929 }
1930 }
1931 }
1932 return m;
1933 }
1934
SBM_read_in_file(SparseBlockStructuredMatrix * const m,FILE * file)1935 void SBM_read_in_file(SparseBlockStructuredMatrix* const m, FILE *file)
1936 {
1937 if(! m)
1938 {
1939 fprintf(stderr, "Numerics, SparseBlockStructuredMatrix SBM_read_in_file failed, NULL input.\n");
1940 exit(EXIT_FAILURE);
1941 }
1942 assert(file);
1943 CHECK_IO(fscanf(file, "%d", &(m->nbblocks)));
1944
1945 if(m->nbblocks == 0) return;
1946
1947 CHECK_IO(fscanf(file, "%d", &(m->blocknumber0)));
1948 CHECK_IO(fscanf(file, "%d", &(m->blocknumber1)));
1949
1950 for(unsigned int i = 0; i < m->blocknumber0; i++)
1951 {
1952 CHECK_IO(fscanf(file, "%d", &(m->blocksize0[i])));
1953 }
1954 for(unsigned int i = 0; i < m->blocknumber1; i++)
1955 {
1956 CHECK_IO(fscanf(file, "%d", &(m->blocksize1[i])));
1957 }
1958
1959 unsigned int filled1 = 0, filled2 = 0;
1960 CHECK_IO(fscanf(file, "%d", &(filled1)));
1961 m->filled1 = filled1;
1962 CHECK_IO(fscanf(file, "%d", &(filled2)));
1963 m->filled2 = filled2;
1964
1965 unsigned int index1_dataCurrent = 0;
1966 for(unsigned int i = 0; i < m->filled1; i++)
1967 {
1968 CHECK_IO(fscanf(file, "%d", &(index1_dataCurrent)));
1969 m->index1_data[i] = index1_dataCurrent;
1970 }
1971 unsigned int index2_dataCurrent = 0;
1972 for(unsigned int i = 0; i < m->filled2; i++)
1973 {
1974 CHECK_IO(fscanf(file, "%d", &(index2_dataCurrent)));
1975 m->index2_data[i] = index2_dataCurrent;
1976 }
1977 unsigned int currentRowNumber ;
1978 size_t colNumber;
1979 unsigned int nbRows, nbColumns;
1980 unsigned int blockk;
1981 for(currentRowNumber = 0 ; currentRowNumber < m->filled1 - 1; ++currentRowNumber)
1982 {
1983 for(size_t blockNum = m->index1_data[currentRowNumber];
1984 blockNum < m->index1_data[currentRowNumber + 1]; ++blockNum)
1985 {
1986 assert(blockNum < m->filled2);
1987 colNumber = m->index2_data[blockNum];
1988 /* Get dim. of the current block */
1989 nbRows = m->blocksize0[currentRowNumber];
1990
1991 if(currentRowNumber != 0)
1992 nbRows -= m->blocksize0[currentRowNumber - 1];
1993
1994 nbColumns = m->blocksize1[colNumber];
1995 if(colNumber != 0)
1996 nbColumns -= m->blocksize1[colNumber - 1];
1997 //fprintf(file,"block[%i] of size %dX%d\n", blockNum, nbRows,nbColumns);
1998
1999 CHECK_IO(fscanf(file, "%d", &(blockk)));
2000 if(blockk != blockNum)
2001 {
2002 printf("Numerics, SparseBlockStructuredMatrix SBM_read_in_file failed, problem in block numbering. \n");
2003 }
2004
2005 for(unsigned int i = 0; i < nbRows * nbColumns; i++)
2006 {
2007 CHECK_IO(fscanf(file, "%32le\n", &(m->block[blockNum][i])));
2008 }
2009
2010 }
2011 }
2012 }
SBM_read_in_filename(SparseBlockStructuredMatrix * const m,const char * filename)2013 void SBM_read_in_filename(SparseBlockStructuredMatrix* const m, const char *filename)
2014 {
2015
2016 }
SBM_clear_pred(SparseBlockStructuredMatrixPred * blmatpred)2017 void SBM_clear_pred(SparseBlockStructuredMatrixPred *blmatpred)
2018 {
2019
2020 for(int i = 0 ; i < blmatpred->nbbldiag ; i++)
2021 {
2022 free(blmatpred->indic[i]);
2023 free(blmatpred->indicop[i]);
2024 free(blmatpred->submatlcp[i]);
2025 free(blmatpred->submatlcpop[i]);
2026 free(blmatpred->ipiv[i]);
2027 free(blmatpred->subq[i]);
2028 free(blmatpred->bufz[i]);
2029 free(blmatpred->newz[i]);
2030 free(blmatpred->workspace[i]);
2031 }
2032 free(blmatpred->indic);
2033 free(blmatpred->indicop);
2034 free(blmatpred->submatlcp);
2035 free(blmatpred->submatlcpop);
2036 free(blmatpred->ipiv);
2037 free(blmatpred->sizesublcp);
2038 free(blmatpred->sizesublcpop);
2039 free(blmatpred->subq);
2040 free(blmatpred->bufz);
2041 free(blmatpred->newz);
2042 free(blmatpred->workspace);
2043
2044 }
2045
2046
SBM_diagonal_block_indices(SparseBlockStructuredMatrix * const M)2047 unsigned int * SBM_diagonal_block_indices(SparseBlockStructuredMatrix* const M)
2048 {
2049 assert(M);
2050 if(M->diagonal_blocks) return M->diagonal_blocks;
2051
2052 unsigned int * diagonal_blocks = (unsigned int*) malloc(M->blocksize0[M->blocknumber0-1]* sizeof(unsigned int));
2053 for(size_t currentRowNumber = 0 ; currentRowNumber < M->filled1 - 1; ++currentRowNumber)
2054 {
2055 for(size_t blockNum = M->index1_data[currentRowNumber];
2056 blockNum < M->index1_data[currentRowNumber + 1]; ++blockNum)
2057 {
2058 if(M->index2_data[blockNum] == currentRowNumber)
2059 {
2060 diagonal_blocks[currentRowNumber] = blockNum;
2061 break;
2062 }
2063 /* Warning: no check if the diagonal block does not exists */
2064 }
2065 }
2066 M->diagonal_blocks = diagonal_blocks;
2067 return diagonal_blocks;
2068 }
2069
2070
SBM_diagonal_block_index(SparseBlockStructuredMatrix * const M,unsigned int row)2071 unsigned int SBM_diagonal_block_index(SparseBlockStructuredMatrix* const M, unsigned int row)
2072 {
2073
2074 unsigned int * diagonal_blocks = SBM_diagonal_block_indices(M);
2075 return diagonal_blocks[row];
2076
2077 /* /\* Look for the first block of row number num *\/ */
2078 /* unsigned int pos; */
2079 /* size_t firstBlockOfRow = M->index1_data[row]; */
2080
2081 /* /\* Look at the diagonal block *\/ */
2082 /* for (pos = (unsigned int)firstBlockOfRow; M->index2_data[pos] != row; ++pos, assert(pos < M->filled2)); */
2083
2084 /* return pos; */
2085 }
2086
SBM_entry(SparseBlockStructuredMatrix * M,unsigned int row,unsigned int col,double val)2087 int SBM_entry(SparseBlockStructuredMatrix* M, unsigned int row, unsigned int col, double val)
2088 {
2089 DEBUG_BEGIN("SBM_entry(...)\n");
2090 DEBUG_PRINTF("row= %i, col =% i, val =%e\n", (int)row, (int)col, val);
2091
2092
2093 /* find the correct row number of blocks */
2094 unsigned int rowNumber =0;
2095 while(M->blocksize0[rowNumber] <= row)
2096 {
2097 rowNumber++;
2098 if(rowNumber >= M->blocknumber0)
2099 {
2100 numerics_warning("SBM_entry", "The row number exceeds the size of the matrix");
2101 DEBUG_END("SBM_entry(...)\n");
2102 return 0;
2103 }
2104 }
2105 DEBUG_PRINTF("rowNumber= %i\n", (int)rowNumber);
2106
2107
2108 unsigned int colNumber =0;
2109
2110 while(M->blocksize1[colNumber] <= col)
2111 {
2112 colNumber++;
2113 if(colNumber >= M->blocknumber1)
2114 {
2115 numerics_warning("SBM_entry", "The col number exceeds the size of the matrix");
2116 DEBUG_END("SBM_entry(...)\n");
2117 return 0;
2118 }
2119 }
2120 DEBUG_PRINTF("colNumber= %i\n", colNumber);
2121
2122 for(size_t blockNum = M->index1_data[rowNumber];
2123 blockNum < M->index1_data[rowNumber + 1]; ++blockNum)
2124 {
2125 assert(blockNum < M->filled2);
2126 if(colNumber == M->index2_data[blockNum])
2127 {
2128 DEBUG_PRINTF("blockNum = %zu\n", blockNum);
2129 /* Get dim. of the current block */
2130 size_t nbRows = M->blocksize0[rowNumber];
2131 size_t row_pos = row;
2132 if(rowNumber != 0)
2133 {
2134 nbRows -= M->blocksize0[rowNumber - 1];
2135 row_pos -= M->blocksize0[rowNumber - 1];
2136 }
2137 assert(nbRows);
2138 assert(row_pos < nbRows);
2139
2140 size_t nbColumns = M->blocksize1[colNumber];
2141 size_t col_pos=col;
2142 if(colNumber != 0)
2143 {
2144 nbColumns -= M->blocksize1[colNumber - 1];
2145 col_pos -= M->blocksize1[colNumber - 1];
2146 }
2147 assert(nbColumns);
2148 assert(col_pos < nbColumns);
2149
2150 DEBUG_PRINTF("row_pos = %zu, col_pos = %zu \n", row_pos, col_pos);
2151
2152 M->block[blockNum][row_pos + col_pos * nbRows] = val;
2153 DEBUG_END("SBM_entry(...)\n");
2154 SBM_inc_version(M);
2155 return 1;
2156 }
2157 }
2158 numerics_warning("SBM_entry", "no existing block for inserting entry\n");
2159 DEBUG_END("SBM_entry(...)\n");
2160 return 0;
2161 }
2162
SBM_get_value(const SparseBlockStructuredMatrix * const M,unsigned int row,unsigned int col)2163 double SBM_get_value(const SparseBlockStructuredMatrix* const M, unsigned int row, unsigned int col)
2164 {
2165 /* Search the row of blocks and the column of blocks */
2166 unsigned int rowblocknumber = M->blocknumber0;
2167 unsigned int colblocknumber = M->blocknumber1;
2168 int rowpos = -1;
2169 int colpos = -1;
2170 int blockNumber = -1;
2171 size_t colNumber;
2172 assert(row < M->blocksize0[M->blocknumber0 - 1]);
2173 assert(col < M->blocksize1[M->blocknumber1 - 1]);
2174
2175
2176 for(unsigned int i = 0; i < M->blocknumber0; i++)
2177 {
2178 if((row < M->blocksize0[i]))
2179 {
2180 rowblocknumber = i;
2181 assert(rowblocknumber < M->blocknumber0);
2182 rowpos = row;
2183 if(i != 0)
2184 rowpos -= M->blocksize0[i - 1];
2185 break;
2186 }
2187
2188 }
2189 for(unsigned int j = 0; j < M->blocknumber1; j++)
2190 {
2191 if((col < M->blocksize1[j]))
2192 {
2193 colblocknumber = j;
2194 assert(colblocknumber < M->blocknumber1);
2195 colpos = col;
2196 if(j != 0)
2197 colpos -= M->blocksize1[j - 1];
2198 blockNumber = -1;
2199 for(size_t blockNum = M->index1_data[rowblocknumber];
2200 blockNum < M->index1_data[rowblocknumber + 1]; ++blockNum)
2201 {
2202 assert(blockNum < M->filled2);
2203
2204 colNumber = M->index2_data[blockNum];
2205 if(colblocknumber == colNumber)
2206 {
2207 blockNumber = (int)blockNum;
2208 break;
2209 }
2210
2211 }
2212 break;
2213 }
2214 }
2215
2216 if(blockNumber == -1) return 0.0;
2217
2218 /* Number of rows/columns of the current block */
2219 int nbRows, nbColumns;
2220 nbRows = M->blocksize0[rowblocknumber];
2221 if(rowblocknumber != 0)
2222 nbRows -= M->blocksize0[rowblocknumber - 1];
2223 nbColumns = M->blocksize1[colblocknumber];
2224 if(colblocknumber != 0)
2225 nbColumns -= M->blocksize1[colblocknumber - 1];
2226 assert(rowpos < nbRows);
2227 assert(colpos < nbColumns);
2228
2229 return M->block[blockNumber][rowpos + colpos * nbRows];
2230
2231
2232
2233 }
2234
SBM_copy(const SparseBlockStructuredMatrix * const A,SparseBlockStructuredMatrix * B,unsigned int copyBlock)2235 int SBM_copy(const SparseBlockStructuredMatrix* const A, SparseBlockStructuredMatrix* B, unsigned int copyBlock)
2236 {
2237 assert(A);
2238 assert(B);
2239
2240 int need_blocks = 0;
2241
2242 if(B->nbblocks < A->nbblocks)
2243 {
2244 need_blocks = 1;
2245 for(unsigned i=0; i<B->nbblocks; ++i)
2246 {
2247 free(B->block [i]);
2248 B->block [i] = NULL;
2249 }
2250 B->block = (double **) realloc(B->block, A->nbblocks * sizeof(double *));
2251 }
2252 B->nbblocks = A->nbblocks;
2253
2254 if(B->blocknumber0 < A->blocknumber0)
2255 {
2256 B->blocksize0 = (unsigned int*) realloc(B->blocksize0, A->blocknumber0 * sizeof(unsigned int));
2257 }
2258 B->blocknumber0 = A->blocknumber0;
2259
2260 if(B->blocknumber1 < A->blocknumber1)
2261 {
2262 B->blocksize1 = (unsigned int*) realloc(B->blocksize1, A->blocknumber1 * sizeof(unsigned int));
2263 }
2264 B->blocknumber1 = A->blocknumber1;
2265
2266 if(B->filled1 < A->filled1)
2267 {
2268 B->index1_data = (size_t*) realloc(B->index1_data, A->filled1 * sizeof(size_t));
2269 }
2270 B->filled1 = A->filled1;
2271
2272 if(B->filled2 < A->filled2)
2273 {
2274 B->index2_data = (size_t*) realloc(B->index2_data, A->filled2 * sizeof(size_t));
2275 }
2276 B->filled2 = A->filled2;
2277
2278
2279 memcpy(B->blocksize0, A->blocksize0, A->blocknumber0 * sizeof(unsigned int));
2280 memcpy(B->blocksize1, A->blocksize1, A->blocknumber1 * sizeof(unsigned int));
2281 memcpy(B->index1_data, A->index1_data, A->filled1 * sizeof(size_t));
2282 memcpy(B->index2_data, A->index2_data, A->filled2 * sizeof(size_t));
2283
2284 if(copyBlock)
2285 {
2286 unsigned int currentRowNumber ;
2287 size_t colNumber;
2288 unsigned int nbRows, nbColumns;
2289 for(currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
2290 {
2291 for(size_t blockNum = A->index1_data[currentRowNumber];
2292 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
2293 {
2294 assert(blockNum < A->filled2);
2295 colNumber = A->index2_data[blockNum];
2296 /* Get dim. of the current block */
2297 nbRows = A->blocksize0[currentRowNumber];
2298 if(currentRowNumber != 0)
2299 nbRows -= A->blocksize0[currentRowNumber - 1];
2300 nbColumns = A->blocksize1[colNumber];
2301
2302 if(colNumber != 0)
2303 nbColumns -= A->blocksize1[colNumber - 1];
2304
2305 if(need_blocks)
2306 {
2307 B->block[blockNum] = (double*)malloc(nbRows * nbColumns * sizeof(double));
2308 }
2309
2310 for(unsigned int i = 0; i < nbRows * nbColumns; i++)
2311 {
2312 B->block[blockNum] [i] = A->block[blockNum] [i] ;
2313 }
2314 }
2315 }
2316 }
2317 else
2318 {
2319 for(unsigned int n = 0; n < B->nbblocks; n++)
2320 B->block[n] = A->block[n];
2321 }
2322
2323 if(B->diagonal_blocks)
2324 {
2325 free(B->diagonal_blocks);
2326 B->diagonal_blocks=NULL;
2327 }
2328
2329 if(A->diagonal_blocks)
2330 {
2331 B->diagonal_blocks = (unsigned int*) malloc(A->blocksize0[A->blocknumber0-1]* sizeof(unsigned int));
2332 memcpy(B->diagonal_blocks,A->diagonal_blocks, A->blocksize0[A->blocknumber0-1]* sizeof(unsigned int));
2333 }
2334
2335
2336
2337
2338 return 0;
2339 }
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
SBM_transpose(const SparseBlockStructuredMatrix * const A,SparseBlockStructuredMatrix * B)2351 int SBM_transpose(const SparseBlockStructuredMatrix* const A, SparseBlockStructuredMatrix* B)
2352 {
2353 assert(A);
2354 assert(B);
2355 B->nbblocks = A->nbblocks;
2356 B->blocknumber0 = A->blocknumber1;
2357 B->blocknumber1 = A->blocknumber0;
2358 B->blocksize0 = (unsigned int*)malloc(B->blocknumber0 * sizeof(unsigned int));
2359 for(unsigned int i = 0; i < B->blocknumber0; i++) B->blocksize0[i] = A->blocksize1[i];
2360 B->blocksize1 = (unsigned int*)malloc(B->blocknumber1 * sizeof(unsigned int));
2361 for(unsigned int i = 0; i < B->blocknumber1; i++) B->blocksize1[i] = A->blocksize0[i];
2362
2363 B->filled1 = A->blocknumber1 + 1;;
2364 B->filled2 = A->filled2;
2365 B->index1_data = (size_t *)malloc(B->filled1 * sizeof(size_t));
2366 B->index2_data = (size_t *)malloc(B->filled2 * sizeof(size_t));
2367
2368
2369 unsigned int currentRowNumberofA;
2370 unsigned int currentColNumberofA ;
2371 size_t colNumberofA;
2372
2373 size_t * blockMap = (size_t *)malloc(B->filled2 * sizeof(size_t));
2374
2375
2376
2377 int blockNumB = -1;
2378 B->index1_data[0] = 0;
2379 for(currentColNumberofA = 0 ; currentColNumberofA < A->blocknumber1; ++currentColNumberofA)
2380 {
2381 assert(currentColNumberofA + 1 < B->filled1);
2382 B->index1_data[currentColNumberofA + 1] = B->index1_data[currentColNumberofA];
2383 for(currentRowNumberofA = 0 ; currentRowNumberofA < A->blocknumber0; ++currentRowNumberofA)
2384 {
2385 for(size_t blockNum = A->index1_data[currentRowNumberofA];
2386 blockNum < A->index1_data[currentRowNumberofA + 1]; ++blockNum)
2387 {
2388 assert(blockNum < B->filled2);
2389 colNumberofA = A->index2_data[blockNum];
2390 if(colNumberofA == currentColNumberofA)
2391 {
2392 blockNumB++;
2393 assert(blockNumB < (int) B->nbblocks && blockNumB >= 0);
2394 B->index1_data[currentColNumberofA + 1]++;
2395 B->index2_data[blockNumB] = currentRowNumberofA;
2396 blockMap[blockNumB] = blockNum;
2397
2398 }
2399 }
2400 }
2401 }
2402
2403 #ifdef VERBOSE_DEBUG
2404 printf("----------------- blockMap ---------------\n");
2405 for(int i = 0; i < B->filled2; i++)
2406 {
2407 printf("blockMap[%i] = %i\n", i, (int)blockMap[i]);
2408 }
2409 printf("----------------- blockMap ---------------\n");
2410 #endif
2411
2412
2413
2414 B->block = (double **)malloc(B->nbblocks * sizeof(double*));
2415 unsigned int currentRowNumber ;
2416 size_t colNumber;
2417 unsigned int nbRows, nbColumns;
2418 for(currentRowNumber = 0 ; currentRowNumber < B->filled1 - 1; ++currentRowNumber)
2419 {
2420 for(size_t blockNum = B->index1_data[currentRowNumber];
2421 blockNum < B->index1_data[currentRowNumber + 1]; ++blockNum)
2422 {
2423 assert(blockNum < B->filled2);
2424 colNumber = B->index2_data[blockNum];
2425 /* Get dim. of the current block */
2426 nbRows = B->blocksize0[currentRowNumber];
2427 if(currentRowNumber != 0)
2428 nbRows -= B->blocksize0[currentRowNumber - 1];
2429 nbColumns = B->blocksize1[colNumber];
2430 if(colNumber != 0)
2431 nbColumns -= B->blocksize1[colNumber - 1];
2432 int lengthblock = nbRows * nbColumns;
2433 B->block[blockNum] = (double*)malloc(lengthblock * sizeof(double));
2434
2435
2436
2437
2438 for(unsigned int i = 0; i < nbRows; i++)
2439 {
2440 for(unsigned int j = 0; j < nbColumns; j++)
2441 {
2442 assert(i + j * nbRows < nbRows * nbColumns);
2443 assert(j + i * nbColumns < nbRows * nbColumns);
2444
2445 B->block[blockNum] [i + j * nbRows ] = A->block[blockMap[blockNum]] [j + i * nbColumns] ;
2446 }
2447 }
2448
2449 }
2450 }
2451 free(blockMap);
2452
2453
2454 return 0;
2455 }
SBM_inverse_diagonal_block_matrix_in_place(const SparseBlockStructuredMatrix * M,int * ipiv)2456 int SBM_inverse_diagonal_block_matrix_in_place(const SparseBlockStructuredMatrix* M, int* ipiv)
2457 {
2458 for(unsigned int i = 0; i < M->filled1 - 1; i++)
2459 {
2460 size_t numberofblockperrow = M->index1_data[i + 1] - M->index1_data[i];
2461 if(numberofblockperrow != 1)
2462 {
2463 fprintf(stderr, "SparseBlockMatrix : SBM_inverse_diagonal_block_matrix: Not a diagonal block matrix\n");
2464 exit(EXIT_FAILURE);
2465 }
2466 }
2467 for(unsigned int i = 0; i < M->filled2; i++)
2468 {
2469 if(M->index2_data[i] != i)
2470 {
2471 fprintf(stderr, "SparseBlockMatrix : SBM_inverse_diagonal_block_matrix: Not a diagonal block matrix\n");
2472 exit(EXIT_FAILURE);
2473 }
2474 }
2475
2476 unsigned int currentRowNumber ;
2477 size_t colNumber;
2478 unsigned int nbRows, nbColumns;
2479 lapack_int infoDGETRF = 0;
2480 lapack_int infoDGETRI = 0;
2481 int info = 0;
2482
2483 lapack_int* lapack_ipiv = (lapack_int *) ipiv;
2484
2485
2486 for(currentRowNumber = 0 ; currentRowNumber < M->filled1 - 1; ++currentRowNumber)
2487 {
2488 for(size_t blockNum = M->index1_data[currentRowNumber];
2489 blockNum < M->index1_data[currentRowNumber + 1]; ++blockNum)
2490 {
2491 assert(blockNum < M->filled2);
2492 colNumber = M->index2_data[blockNum];
2493 /* Get dim. of the current block */
2494 nbRows = M->blocksize0[currentRowNumber];
2495 if(currentRowNumber != 0)
2496 nbRows -= M->blocksize0[currentRowNumber - 1];
2497 nbColumns = M->blocksize1[colNumber];
2498 if(colNumber != 0)
2499 nbColumns -= M->blocksize1[colNumber - 1];
2500
2501 assert(nbRows == nbColumns);
2502
2503 DGETRF(nbRows, nbColumns, M->block[blockNum], nbRows, lapack_ipiv, &infoDGETRF);
2504 assert(!infoDGETRF);
2505
2506 DGETRI(nbRows, M->block[blockNum], nbRows, lapack_ipiv, &infoDGETRI);
2507 assert(!infoDGETRI);
2508
2509 }
2510 }
2511 if((!infoDGETRF) || (!infoDGETRI)) info = 0;
2512
2513 return info;
2514
2515
2516 }
2517
SBM_to_dense(const SparseBlockStructuredMatrix * const A,double * denseMat)2518 void SBM_to_dense(const SparseBlockStructuredMatrix* const A, double *denseMat)
2519 {
2520 assert(A);
2521 assert(A->blocksize0);
2522 assert(A->blocksize1);
2523 int n = A->blocksize0[A->blocknumber0 - 1];
2524 int m = A->blocksize1[A->blocknumber1 - 1];
2525
2526 /* denseMat = (double*)malloc(n*m*sizeof(double)); */
2527 for(int i = 0; i < n ; i++)
2528 {
2529 for(int j = 0; j < m; j++)
2530 {
2531 denseMat[i + j * n] = SBM_get_value(A, i, j);
2532 }
2533 }
2534 }
2535
SBM_to_sparse_init_memory(const SparseBlockStructuredMatrix * const A,CSparseMatrix * sparseMat)2536 int SBM_to_sparse_init_memory(const SparseBlockStructuredMatrix* const A, CSparseMatrix *sparseMat)
2537 {
2538 DEBUG_BEGIN("SBM_to_sparse_init_memory(...)\n")
2539 assert(A);
2540 assert(A->blocksize0);
2541 assert(A->blocksize1);
2542 int n = A->blocksize0[A->blocknumber0 - 1];
2543 int m = A->blocksize1[A->blocknumber1 - 1];
2544
2545 sparseMat->m = n;
2546 sparseMat->n = m;
2547
2548 sparseMat->nz = -2; /* csr */
2549 sparseMat->nzmax = 0;
2550 sparseMat->p = (CS_INT*)malloc((sparseMat->m + 1) * sizeof(CS_INT));
2551
2552 /* Row (block) position of the current block */
2553 unsigned int currentRowNumber ;
2554 /* Column (block) position of the current block*/
2555 size_t colNumber;
2556 /* Number of rows/columns of the current block */
2557 int nbRows, nbColumns;
2558
2559 for(currentRowNumber = 0 ; currentRowNumber < A->filled1 - 1; ++currentRowNumber)
2560 {
2561 for(size_t blockNum = A->index1_data[currentRowNumber];
2562 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
2563 {
2564 assert(blockNum < A->filled2);
2565 colNumber = A->index2_data[blockNum];
2566
2567 /* Get dim. of the current block */
2568 nbRows = A->blocksize0[currentRowNumber];
2569
2570 if(currentRowNumber != 0)
2571 nbRows -= A->blocksize0[currentRowNumber - 1];
2572 assert((nbRows >= 0));
2573 nbColumns = A->blocksize1[colNumber];
2574 if(colNumber != 0)
2575 nbColumns -= A->blocksize1[colNumber - 1];
2576 assert((nbColumns >= 0));
2577
2578 sparseMat->nzmax += nbColumns * nbRows;
2579 }
2580 }
2581 sparseMat->i = (CS_INT*)malloc((sparseMat->nzmax) * sizeof(CS_INT));
2582 sparseMat->x = (double*)malloc((sparseMat->nzmax) * sizeof(double));
2583 DEBUG_END("SBM_to_sparse_init_memory(...)\n")
2584 return 0;
2585 }
2586
sparseMatrixBegin(const CSparseMatrix * const sparseMat)2587 sparse_matrix_iterator sparseMatrixBegin(const CSparseMatrix* const sparseMat)
2588 {
2589 if(sparseMat->nz >= 0)
2590 {
2591 sparse_matrix_iterator i = { 0, -1, 0, 0, NAN, sparseMat};
2592 return i;
2593 }
2594 else if(sparseMat->nz == -1)
2595 {
2596 sparse_matrix_iterator i = { 0, sparseMat->p[0], 0, 0, NAN, sparseMat };
2597 return i;
2598 }
2599 else if(sparseMat->nz == -2)
2600 {
2601 sparse_matrix_iterator i = { 0, sparseMat->p[0], 0, 0, NAN, sparseMat };
2602 return i;
2603 }
2604
2605 assert(0);
2606 sparse_matrix_iterator z = { 0, 0, 0, 0, NAN, sparseMat };
2607 return z;
2608 }
2609
sparseMatrixNext(sparse_matrix_iterator * it)2610 int sparseMatrixNext(sparse_matrix_iterator* it)
2611 {
2612 if(it->mat->nz >= 0) /* triplet */
2613 {
2614 assert(it->counter2 == -1);
2615
2616 if(it->counter1 < it->mat->nz)
2617 {
2618 assert(it->counter1 >= 0);
2619
2620 it->first = it->mat->i[it->counter1];
2621 it->second = it->mat->p[it->counter1];
2622 it->third = it->mat->x[it->counter1];
2623
2624 it->counter1++;
2625
2626 return 1;
2627 }
2628 else
2629 {
2630 return 0; /* stop */
2631 }
2632 }
2633 else if(it->mat->nz == -1) /* csc */
2634 {
2635 if(it->counter1 < it->mat->n)
2636 {
2637 if(it->counter2 < it->mat->p[it->counter1 + 1])
2638 {
2639 /* next line */
2640 assert(it->counter2 >= 0);
2641 assert(it->counter2 < it->mat->nzmax);
2642
2643 it->first = it->mat->i[it->counter2];
2644 it->second = it->counter1;
2645 it->third = it->mat->x[it->counter2];
2646 it->counter2++;
2647 return 1;
2648 }
2649 else
2650 {
2651 /* next column */
2652 it->counter1++;
2653
2654 assert(it->counter1 >= 0);
2655 assert(it->counter1 < (it->mat->n + 1));
2656 it->counter2 = it->mat->p[it->counter1];
2657
2658 DEBUG_PRINTF("it->counter1 = %i, it->counter2 = %i\n", (int)it->counter1, (int)it->counter2);
2659
2660 if(it->counter2 < it->mat->nzmax)
2661 {
2662
2663 assert(it->counter2 >= 0);
2664 assert(it->counter2 < it->mat->nzmax);
2665
2666 it->first = it->mat->i[it->counter2];
2667 it->second = it->counter1;
2668 it->third = it->mat->x[it->counter2];
2669
2670 it->counter2++;
2671 return 1;
2672 }
2673 else
2674 {
2675 return 0; /* stop */
2676 }
2677 }
2678 }
2679 return 0; /* stop */
2680 }
2681 else if(it->mat->nz == -2) /* csr */
2682 {
2683 if(it->counter1 < it->mat->m)
2684 {
2685 if(it->counter2 < it->mat->p[it->counter1 + 1])
2686 {
2687 /* next column */
2688
2689 assert(it->counter2 >= 0);
2690 assert(it->counter2 < it->mat->nzmax);
2691
2692 it->first = it->counter1;
2693 it->second = it->mat->i[it->counter2];
2694 it->third = it->mat->x[it->counter2];
2695 it->counter2++;
2696 return 1;
2697 }
2698 else
2699 {
2700 /* next line */
2701 it->counter1++;
2702
2703 assert(it->counter1 >= 0);
2704 assert(it->counter1 < (it->mat->m + 1));
2705
2706 it->counter2 = it->mat->p[it->counter1];
2707
2708 DEBUG_PRINTF("it->counter1 = %i, it->counter2 = %i\n", (int)it->counter1, (int)it->counter2);
2709
2710 if(it->counter2 < it->mat->nzmax)
2711 {
2712 assert(it->counter2 >= 0);
2713 assert(it->counter2 < it->mat->nzmax);
2714
2715 it->first = it->counter1;
2716 it->second = it->mat->i[it->counter2];
2717 it->third = it->mat->x[it->counter2];
2718
2719 it->counter2++;
2720 return 1;
2721 }
2722 else
2723 {
2724 return 0; /* stop */
2725 }
2726 }
2727 }
2728 else
2729 {
2730 return 0; /* stop */
2731 }
2732 }
2733
2734 assert(0);
2735 return 0;
2736 }
2737
2738
SBCM_null(SparseBlockCoordinateMatrix * MC)2739 void SBCM_null(SparseBlockCoordinateMatrix* MC)
2740 {
2741 MC->blocknumber0 = 0;
2742 MC->blocknumber1 = 0;
2743 MC->nbblocks = 0;
2744
2745 MC->row = NULL;
2746 MC->column = NULL;
2747 MC->blocksize0 = NULL;
2748 MC->blocksize1 = NULL;
2749 MC->block=NULL ;
2750
2751 }
2752
SBCM_new(void)2753 SparseBlockCoordinateMatrix* SBCM_new(void)
2754 {
2755 SparseBlockCoordinateMatrix* MC = (SparseBlockCoordinateMatrix*) malloc(sizeof(SparseBlockCoordinateMatrix));
2756 SBCM_null(MC);
2757 return MC;
2758 }
2759
2760
2761
SBCM_new_3x3(unsigned int m,unsigned int n,unsigned int nbblocks,unsigned int * row,unsigned int * column,double * block)2762 SparseBlockCoordinateMatrix* SBCM_new_3x3(unsigned int m, unsigned int n,
2763 unsigned int nbblocks,
2764 unsigned int *row,
2765 unsigned int *column,
2766 double *block)
2767 {
2768
2769 SparseBlockCoordinateMatrix* MC = SBCM_new();
2770 MC->blocknumber0 = m; /* block row */
2771 MC->blocknumber1 = n; /* block column */
2772 MC->nbblocks = nbblocks;
2773 MC->row = (unsigned int *) malloc(sizeof(unsigned int) * nbblocks);
2774 MC->column = (unsigned int *) malloc(sizeof(unsigned int) * nbblocks);
2775 for(unsigned int i = 0; i < nbblocks; ++i)
2776 {
2777 MC->row[i] = row[i] - 1;
2778 }
2779 for(unsigned int i = 0; i < nbblocks; ++i)
2780 {
2781 MC->column[i] = column[i] - 1;
2782 }
2783 MC->block = (double **) malloc(sizeof(double*)*nbblocks);
2784 MC->blocksize0 = (unsigned int *) malloc(sizeof(unsigned int) * nbblocks);
2785 MC->blocksize1 = (unsigned int *) malloc(sizeof(unsigned int) * nbblocks);
2786 for(unsigned int i = 0; i < nbblocks; ++i)
2787 {
2788 MC->blocksize0[i] = 3 * (i + 1);
2789 MC->blocksize1[i] = 3 * (i + 1);
2790 MC->block[i] = &block[i * 9];
2791 }
2792
2793 return MC;
2794 }
2795
SBCM_free_3x3(SparseBlockCoordinateMatrix * MC)2796 void SBCM_free_3x3(SparseBlockCoordinateMatrix *MC)
2797 {
2798 free(MC->block);
2799 free(MC->blocksize0);
2800 free(MC->blocksize1);
2801 free(MC->row);
2802 free(MC->column);
2803 }
2804
2805 /* quite obvious alg but in case of incomprehension see Nathan Bell coo_tocsr */
2806 /* i.e coo.h file under scipy sparsetools */
SBCM_to_SBM(SparseBlockCoordinateMatrix * MC)2807 SparseBlockStructuredMatrix* SBCM_to_SBM(SparseBlockCoordinateMatrix* MC)
2808 {
2809 SparseBlockStructuredMatrix* M = SBM_new();
2810
2811 M->nbblocks = MC->nbblocks;
2812 M->filled2 = MC->nbblocks;
2813 M->blocknumber0 = MC->blocknumber0;
2814 M->blocknumber1 = MC->blocknumber1;
2815 M->blocksize0 = MC->blocksize0;
2816 M->blocksize1 = MC->blocksize1;
2817
2818 M->index1_data = (size_t *) calloc(M->blocknumber0 + 1, sizeof(size_t));
2819 M->index2_data = (size_t *) malloc(sizeof(size_t) * M->nbblocks);
2820 M->filled1 = 2;
2821 for(unsigned int i = 0; i < M->nbblocks; ++i)
2822 {
2823 M->index1_data[MC->row[i] + 1]++;
2824 M->filled1 = ((MC->row[i] + 2) > M->filled1) ? MC->row[i] + 2 : M->filled1;
2825 }
2826
2827 for(unsigned int i = 1; i < M->blocknumber0 + 1; ++i)
2828 {
2829 M->index1_data[i] += M->index1_data[i - 1];
2830 }
2831
2832 M->block = (double **) malloc(sizeof(double*)*M->nbblocks);
2833
2834 for(unsigned int i = 0; i < M->nbblocks; ++i)
2835 {
2836 unsigned int row = MC->row[i];
2837 size_t dest = M->index1_data[row];
2838
2839 assert(dest < M->nbblocks);
2840 M->index2_data[dest] = MC->column[i];
2841
2842 assert(dest < M->nbblocks);
2843 M->block[dest] = MC->block[i];
2844
2845 assert(row < (M->blocknumber0 + 1));
2846 M->index1_data[row]++;
2847 }
2848
2849 size_t last = 0;
2850 for(unsigned int i = 0; i <= M->blocknumber0; i++)
2851 {
2852 size_t temp = M->index1_data[i];
2853 M->index1_data[i] = last;
2854 last = temp;
2855 }
2856
2857 return M;
2858 }
2859
SBM_free_from_SBCM(SparseBlockStructuredMatrix * M)2860 void SBM_free_from_SBCM(SparseBlockStructuredMatrix* M)
2861 {
2862 free(M->index1_data);
2863 free(M->index2_data);
2864 free(M->block);
2865 free(M);
2866 }
2867
SBM_from_csparse(int blocksize,const CSparseMatrix * const sparseMat,SparseBlockStructuredMatrix * A)2868 int SBM_from_csparse(int blocksize, const CSparseMatrix* const sparseMat, SparseBlockStructuredMatrix* A)
2869 {
2870 DEBUG_PRINT("SBM_from_csparse start\n")
2871 assert(sparseMat);
2872 assert(sparseMat->p);
2873 assert(sparseMat->i);
2874 assert(sparseMat->x);
2875 /* assert(sparseMat->nz == -2); */
2876
2877 assert(sparseMat->m % blocksize == 0);
2878 assert(sparseMat->n % blocksize == 0);
2879
2880 CS_INT bnrow = sparseMat->m / blocksize;
2881 CS_INT bncol = sparseMat->n / blocksize;
2882 DEBUG_PRINTF("SBM_from_csparse. bnrow =%i\n", (int)bnrow);
2883 DEBUG_PRINTF("SBM_from_csparse. bncol =%i\n", (int)bncol);
2884 A->blocknumber0 = (int) bnrow;
2885 A->blocknumber1 = (int) bncol;
2886
2887
2888 A->diagonal_blocks=NULL;
2889
2890 // assert(A->blocksize0 == NULL);
2891 A->blocksize0 = (unsigned int*) malloc(A->blocknumber0 * sizeof(unsigned int));
2892
2893 // assert(A->blocksize1 == NULL);
2894 A->blocksize1 = (unsigned int*) malloc(A->blocknumber1 * sizeof(unsigned int));
2895
2896 for(unsigned int i = 0; i < A->blocknumber0; i++)
2897 {
2898 A->blocksize0[i] = (i + 1) * blocksize;
2899 }
2900
2901 for(unsigned int i = 0; i < A->blocknumber1; i++)
2902 {
2903 A->blocksize1[i] = (i + 1) * blocksize;
2904 }
2905
2906 /* we have to find non empty blocks */
2907
2908 CS_INT blockindexmax = -1;
2909 CS_INT blocklinemax = -1;
2910 int* blockline;
2911 int* blocknum;
2912
2913 /* 1: find blockindexmax (<= bnrow + bncol * bnrow) & blocklinemax */
2914 for(sparse_matrix_iterator it = sparseMatrixBegin(sparseMat);
2915 sparseMatrixNext(&it);)
2916 {
2917
2918 CS_INT row = it.first;
2919 CS_INT col = it.second;
2920
2921 DEBUG_PRINTF("it.first = %i, it.second = %i \n", (int)row, (int)col);
2922 DEBUG_PRINTF("it.third = %g, \n", it.third);
2923
2924 CS_INT brow = row / blocksize;
2925 CS_INT bcol = col / blocksize;
2926 CS_INT blockindex = brow * bncol + bcol;
2927
2928 if((fabs(it.third) > 0.0) && (blockindex > blockindexmax - 1))
2929 {
2930 blockindexmax = blockindex + 1;
2931 };
2932
2933 if((fabs(it.third) > 0.0) && brow > (blocklinemax - 1))
2934 {
2935 blocklinemax = brow + 1;
2936 }
2937 }
2938 DEBUG_PRINTF("SBM_from_csparse. blockindexmax =%i\n", (int)blockindexmax);
2939 DEBUG_PRINTF("SBM_from_csparse. blocklinemax=%i\n", (int)blocklinemax);
2940
2941 // assert(blockindexmax <= bnrow + bncol * bnrow + 1);
2942 // assert(blocklinemax <= bnrow + 1);
2943
2944 /* 2: allocate temporary memory for blocknumbers & blocklines */
2945 blocknum = (int *) malloc(blockindexmax * sizeof(int));
2946 blockline = (int *) malloc(blocklinemax * sizeof(int));
2947 for(int i = 0; i < blockindexmax; i++)
2948 {
2949 blocknum[i] = -1;
2950 }
2951
2952 for(int i = 0; i < blocklinemax; i++)
2953 {
2954 blockline[i] = -1;
2955 }
2956
2957 /* 3: flag non empty blocks & lines */
2958 for(sparse_matrix_iterator it = sparseMatrixBegin(sparseMat);
2959 sparseMatrixNext(&it);)
2960 {
2961 CS_INT row = it.first;
2962 CS_INT col = it.second;
2963
2964 CS_INT brow = row / blocksize;
2965 CS_INT bcol = col / blocksize;
2966
2967 CS_INT blockindex = brow * bncol + bcol;
2968
2969 if(fabs(it.third) > 0.)
2970 {
2971 assert(blockindex < blockindexmax);
2972 blocknum[blockindex] = -2;
2973 assert(brow < blocklinemax);
2974 blockline[brow] = -2;
2975 }
2976 }
2977
2978 /* 4: count non empty blocks */
2979 A->nbblocks = 0;
2980 for(int i = 0; i < blockindexmax; i++)
2981 {
2982 assert(blocknum[i] == -1 || blocknum[i] == -2);
2983
2984 if(blocknum[i] == -2)
2985 {
2986 DEBUG_PRINTF("blocknum[%d] = %d\n", i, (int)A->nbblocks);
2987 blocknum[i] = A->nbblocks++;
2988 }
2989 }
2990
2991 /* 5: allocate memory for contiguous blocks */
2992 assert(A->nbblocks>0);
2993
2994 if(!A->block)
2995 A->block = (double **) malloc(A->nbblocks * sizeof(double *));
2996 for(unsigned int i = 0; i < A->nbblocks; i++)
2997 {
2998 A->block[i] = (double *) malloc(blocksize * blocksize * sizeof(double));
2999
3000 /* fill block with 0 */
3001 for(int birow = 0; birow < blocksize; ++birow)
3002 {
3003 for(int bicol = 0; bicol < blocksize; ++bicol)
3004 {
3005 A->block[i][birow + bicol * blocksize] = 0.;
3006 }
3007 }
3008 }
3009
3010 A->filled2 = (size_t) A->nbblocks; /* one of them should be deprecated! */
3011
3012 /* 6: count non empty lines */
3013 A->filled1 = 1; /* A->filled1 = number of non empty lines + 1 */
3014
3015 for(int i = 0; i < blocklinemax; i++)
3016 {
3017 assert(blockline[i] == -1 || blockline[i] == -2);
3018
3019 if(blockline[i] == -2)
3020 {
3021 A->filled1++;
3022 }
3023 }
3024 DEBUG_PRINTF("A->filled1 =%zu\n",A->filled1);
3025 DEBUG_PRINTF("A->filled2 =%zu\n",A->filled2);
3026 /* 7: allocate memory for index data */
3027
3028 if(!A->index1_data)
3029 A->index1_data = (size_t*) malloc(A->filled1 * sizeof(size_t));
3030
3031 if(!A->index2_data)
3032 A->index2_data = (size_t*) malloc(A->filled2 * sizeof(size_t));
3033
3034
3035 for(size_t i = 0; i < A->filled1; i++)
3036 {
3037 A->index1_data[i] = blockindexmax + 1;
3038 }
3039
3040 for(size_t i = 0; i < A->filled2; i++)
3041 {
3042 A->index2_data[i] = 0;
3043 }
3044
3045 /* 8: fill index1_data & index2_data & copy values in contiguous
3046 * blocks */
3047 for(sparse_matrix_iterator it = sparseMatrixBegin(sparseMat);
3048 sparseMatrixNext(&it);)
3049 {
3050
3051 CS_INT row = it.first;
3052 CS_INT col = it.second;
3053
3054 assert(row < sparseMat->m);
3055 assert(col < sparseMat->n);
3056
3057 CS_INT brow = row / blocksize;
3058 CS_INT bcol = col / blocksize;
3059
3060 assert(brow < bnrow);
3061 assert(bcol < bncol);
3062
3063 CS_INT blockindex = brow * bncol + bcol;
3064
3065 CS_INT birow = row % blocksize; /* block inside row */
3066 CS_INT bicol = col % blocksize; /* block inside column */
3067
3068 A->index1_data[A->filled1 - 1] = A->filled2;
3069
3070 if((blockindex < blockindexmax) && (blocknum[blockindex] >= 0))
3071 {
3072
3073 /* this is an non empty block */
3074
3075 /* index1_data[rowNumber]<= blockNumber */
3076
3077 assert(brow < (CS_INT)A->filled1);
3078 if(A->index1_data[brow] > (size_t)blocknum[blockindex])
3079 {
3080 A->index1_data[brow] = blocknum[blockindex];
3081 }
3082
3083 A->index2_data[blocknum[blockindex]] = bcol;
3084
3085 assert(birow + bicol * blocksize <= blocksize * blocksize);
3086
3087 assert(blockindex < blockindexmax);
3088 assert(blocknum[blockindex] < (CS_INT)A->nbblocks);
3089 assert(blocknum[blockindex] < (CS_INT)A->filled2);
3090
3091 DEBUG_PRINTF("A->block[blocknum[blockindex=%ld]=%d][birow=%ld + bicol=%ld * blocksize=%d] = it.third=%g\n",
3092 blockindex, blocknum[blockindex], birow, bicol, blocksize, it.third);
3093 A->block[blocknum[blockindex]][birow + bicol * blocksize] = it.third;
3094 }
3095 }
3096
3097 /* 9: free temp memory */
3098 free(blocknum);
3099 free(blockline);
3100 DEBUG_PRINT("SBM_from_csparse end\n")
3101
3102 return 0;
3103
3104 }
3105
SBM_to_sparse(const SparseBlockStructuredMatrix * const A,CSparseMatrix * outSparseMat)3106 int SBM_to_sparse(const SparseBlockStructuredMatrix* const A, CSparseMatrix *outSparseMat)
3107 {
3108 DEBUG_BEGIN("SBM_to_sparse(...)\n");
3109 assert(A);
3110 assert(A->blocksize0);
3111 assert(A->blocksize1);
3112
3113 assert(outSparseMat);
3114 assert(outSparseMat->p);
3115 assert(outSparseMat->i);
3116 assert(outSparseMat->x);
3117
3118 /* Row (block) position of the current block */
3119 unsigned int currentRowNumber ;
3120 /* Column (block) position of the current block*/
3121 size_t colNumber;
3122
3123 int nnz = 0;
3124 int isparserowstart;
3125 int isparsecolumnstart;
3126 outSparseMat->p[0] = 0; /* We assume that the first row is non empty */
3127 for(currentRowNumber = 0 ; currentRowNumber < (unsigned int) A->filled1 - 1; ++currentRowNumber)
3128 {
3129
3130 /* Get row dim. of the current block line*/
3131 int nbRows = A->blocksize0[currentRowNumber];
3132
3133 int isparserowend = A->blocksize0[currentRowNumber];
3134
3135 if(currentRowNumber != 0)
3136 {
3137 nbRows -= A->blocksize0[currentRowNumber - 1];
3138 isparserowstart = A->blocksize0[currentRowNumber - 1];
3139 }
3140 else
3141 {
3142 isparserowstart = 0;
3143
3144 }
3145 assert((nbRows >= 0));
3146
3147 DEBUG_PRINTF("isparserowstart = %i\t", isparserowstart);
3148 DEBUG_PRINTF("isparserowend = %i\n", isparserowend);
3149
3150
3151 for(int isparserow = isparserowstart; isparserow < isparserowend ; isparserow++)
3152 {
3153
3154
3155 outSparseMat->p[isparserow + 1] = outSparseMat->p[isparserow];
3156
3157 for(size_t blockNum = A->index1_data[currentRowNumber];
3158 blockNum < A->index1_data[currentRowNumber + 1]; ++blockNum)
3159 {
3160 colNumber = A->index2_data[blockNum];
3161
3162 int nbColumns = A->blocksize1[colNumber];
3163 int isparsecolumnend = A->blocksize1[colNumber];
3164
3165 if(colNumber != 0)
3166 {
3167 nbColumns -= A->blocksize1[colNumber - 1];
3168 isparsecolumnstart = A->blocksize1[colNumber - 1];
3169 }
3170 else
3171 isparsecolumnstart = 0;
3172
3173 assert((nbColumns >= 0));
3174
3175 outSparseMat->p[isparserow + 1] += nbColumns ;
3176
3177
3178 DEBUG_PRINTF("isparsecolumnstart = %i\t", isparsecolumnstart);
3179 DEBUG_PRINTF("isparsecolumnend = %i\n", isparsecolumnend);
3180
3181
3182 for(int isparsecolumn = isparsecolumnstart; isparsecolumn < isparsecolumnend; isparsecolumn++)
3183 {
3184
3185 DEBUG_PRINTF("nnz = %i \t", nnz);
3186 DEBUG_PRINTF("isparserow = %i \t", isparserow);
3187 DEBUG_PRINTF("isparsecolumn = %i \t", isparsecolumn);
3188
3189 outSparseMat->i[nnz] = isparsecolumn;
3190
3191 int rowintheblock = isparserow;
3192 if(currentRowNumber != 0)
3193 {
3194 rowintheblock -= A->blocksize0[currentRowNumber - 1];
3195 }
3196
3197
3198
3199 int colintheblock = isparsecolumn;
3200 if(colNumber != 0)
3201 {
3202 colintheblock -= A->blocksize1[colNumber - 1];
3203 }
3204 assert(rowintheblock < nbRows);
3205 assert(colintheblock < nbColumns);
3206
3207 DEBUG_PRINTF(" rowintheblock= %i \t", rowintheblock);
3208 DEBUG_PRINTF("colintheblock = %i \n", colintheblock);
3209
3210
3211 outSparseMat->x[nnz] = A->block[blockNum][colintheblock * nbRows + rowintheblock];
3212 nnz++;
3213 }
3214 DEBUG_PRINT("\n");
3215
3216 }
3217
3218
3219 }
3220 }
3221 DEBUG_END("SBM_to_sparse(...)\n");
3222 return 0;
3223 }
3224
3225
SBMfree(SparseBlockStructuredMatrix * A,unsigned int level)3226 void SBMfree(SparseBlockStructuredMatrix* A, unsigned int level)
3227 {
3228
3229 if(level & NUMERICS_SBM_FREE_BLOCK)
3230 {
3231 for(unsigned int i = 0; i < A->nbblocks; i++)
3232 free(A->block[i]);
3233 }
3234 free(A->block);
3235 free(A->blocksize0);
3236 free(A->blocksize1);
3237 free(A->index1_data);
3238 free(A->index2_data);
3239 if(level & NUMERICS_SBM_FREE_SBM)
3240 {
3241 printf("val1 : %d", NUMERICS_SBM_FREE_SBM);
3242 printf("val2 : %d", (int)level);
3243 free(A);
3244 }
3245 }
3246
3247 //#define SBM_DEBUG_SBM_row_to_dense
SBM_row_to_dense(const SparseBlockStructuredMatrix * const A,int row,double * denseMat,int rowPos,int rowNb)3248 void SBM_row_to_dense(const SparseBlockStructuredMatrix* const A, int row, double *denseMat, int rowPos, int rowNb)
3249 {
3250 assert(A);
3251 int BlockRowNb = 0;
3252 int ColNb = A->blocksize1[A->blocknumber1 - 1];
3253 if(row)
3254 BlockRowNb = A->blocksize0[row] - A->blocksize0[row - 1];
3255 else
3256 BlockRowNb = A->blocksize0[row];
3257 #ifdef SBM_DEBUG_SBM_row_to_dense
3258 printf("SBM_row_to_dense : copi block row %i, containing %i row and %i col.\n", row, BlockRowNb, ColNb);
3259 #endif
3260
3261 //zero memory
3262 for(int numRow = rowPos; numRow < rowPos + BlockRowNb; numRow++)
3263 for(int numCol = 0; numCol < ColNb; numCol++)
3264 denseMat[numRow + numCol * rowNb] = 0.0;
3265
3266 //index1_data[rowNumber]<= blockNumber <index1_data[rowNumber+1]
3267 for(size_t numBlock = A->index1_data[row]; numBlock < A->index1_data[row + 1]; numBlock++)
3268 {
3269 double * beginBlock = A->block[numBlock];
3270 size_t colNumber = A->index2_data[numBlock];
3271 int indexColBegin = 0;
3272 if(colNumber)
3273 indexColBegin = A->blocksize1[colNumber - 1];
3274 int BlockColNb = A->blocksize1[colNumber] - indexColBegin;
3275 for(int i = 0; i < BlockRowNb; i++)
3276 {
3277 for(int j = 0; j < BlockColNb; j++)
3278 {
3279 denseMat[rowPos + i + (indexColBegin + j)*rowNb] = beginBlock[i + j * BlockRowNb];
3280 }
3281 }
3282 }
3283 #ifdef SBM_DEBUG_SBM_row_to_dense
3284 printf("SBM_row_to_dense : res in file SBM_row_to_dense.txt.");
3285 FILE * titi = fopen("SBM_row_to_dense.txt", "w");
3286 SBM_write_in_fileForScilab(A, titi);
3287 fprintf(titi, "\n//Dense matrix of row block %i:\n", row);
3288 fprintf(titi, "denseRow = [ \t");
3289 for(int i = 0; i < BlockRowNb; i++)
3290 {
3291 fprintf(titi, "[");
3292 for(int j = 0; j < ColNb; j++)
3293 {
3294 fprintf(titi, "%32.24e\t ", denseMat[rowPos + i + j * rowNb]);
3295 }
3296 fprintf(titi, "];\n");
3297 }
3298 fprintf(titi, "];\n");
3299 fclose(titi);
3300 #endif
3301 }
3302 //#define SBM_DEBUG_SBM_ROW_PERM
SBM_row_permutation(unsigned int * rowIndex,SparseBlockStructuredMatrix * A,SparseBlockStructuredMatrix * C)3303 void SBM_row_permutation(unsigned int *rowIndex, SparseBlockStructuredMatrix* A, SparseBlockStructuredMatrix* C)
3304 {
3305 #ifdef SBM_DEBUG_SBM_ROW_PERM
3306 FILE * titi = fopen("SBM_row_permutation_input.txt", "w");
3307 SBM_write_in_fileForScilab(A, titi);
3308 fclose(titi);
3309 #endif
3310 int nbRow = A->blocknumber0;
3311 int nbCol = A->blocknumber1;
3312 C->nbblocks = A->nbblocks;
3313 C->block = (double**)malloc(A->nbblocks * sizeof(double*));
3314 C->blocknumber0 = A->blocknumber0;
3315 C->blocknumber1 = A->blocknumber1;
3316 C->blocksize0 = (unsigned int*)malloc(nbRow * sizeof(unsigned int));
3317 C->blocksize1 = (unsigned int*)malloc(nbCol * sizeof(unsigned int));
3318 C->filled1 = A->filled1;
3319 C->filled2 = A->filled2;
3320 C->index1_data = (size_t*)malloc(C->filled1 * sizeof(size_t));
3321 C->index2_data = (size_t*)malloc(C->filled2 * sizeof(size_t));
3322 /*Row permutation ==> same col size*/
3323 for(int i = 0; i < nbCol; i++)
3324 {
3325 C->blocksize1[i] = A->blocksize1[i];
3326 }
3327 int curNbBlockC = 0;
3328 C->index1_data[0] = 0;
3329 for(int rowC = 0; rowC < nbCol; rowC++)
3330 {
3331 int rowA = rowIndex[rowC];
3332 /*C->blocksize0[rowC+1]=C->blocksize0[rowC]+ number of row of the current block*/
3333 int nbRowInBlock = A->blocksize0[rowA];
3334 if(rowA)
3335 nbRowInBlock -= A->blocksize0[rowA - 1];
3336 #ifdef SBM_DEBUG_SBM_ROW_PERM
3337 printf("SBM_row_permutation rowA=%i, rowC=%i\n", rowA, rowC);
3338 #endif
3339 if(rowC)
3340 C->blocksize0[rowC] = C->blocksize0[rowC - 1] + nbRowInBlock;
3341 else
3342 C->blocksize0[rowC] = nbRowInBlock;
3343 size_t NbBlockCurRow = A->index1_data[rowA + 1] - A->index1_data[rowA];
3344
3345 C->index1_data[rowC + 1] = C->index1_data[rowC] + NbBlockCurRow;
3346 for(size_t numBlockInRowA = A->index1_data[rowA]; numBlockInRowA < A->index1_data[rowA + 1]; numBlockInRowA++)
3347 {
3348 C->index2_data[curNbBlockC] = A->index2_data[numBlockInRowA];
3349 C->block[curNbBlockC] = A->block[numBlockInRowA];
3350 curNbBlockC++;
3351 }
3352 }
3353 #ifdef SBM_DEBUG_SBM_ROW_PERM
3354 titi = fopen("SBM_row_permutation_output.txt", "w");
3355 SBM_write_in_fileForScilab(C, titi);
3356 fclose(titi);
3357 #endif
3358 }
3359 //#define SBM_DEBUG_SBM_COL_PERM
SBM_column_permutation(unsigned int * colIndex,SparseBlockStructuredMatrix * A,SparseBlockStructuredMatrix * C)3360 void SBM_column_permutation(unsigned int *colIndex, SparseBlockStructuredMatrix* A, SparseBlockStructuredMatrix* C)
3361 {
3362 #ifdef SBM_DEBUG_SBM_COL_PERM
3363 FILE * titi = fopen("SBM_column_permutation_input.txt", "w");
3364 SBM_write_in_fileForScilab(A, titi);
3365 fclose(titi);
3366 #endif
3367 SBM_copy(A, C, 0);
3368 for(unsigned int n = 0; n < C->nbblocks; n++)
3369 {
3370 C->index2_data[n] = colIndex[C->index2_data[n]];
3371 }
3372 //int curColnb=0;
3373 int nbBlockCol = A->blocknumber1;
3374 for(int numCol = 0; numCol < nbBlockCol; numCol++)
3375 {
3376 #ifdef SBM_DEBUG_SBM_COL_PERM
3377 printf("SBM_column_permutation colA=%i, colC=%i\n", numCol, (int)colIndex[numCol]);
3378 #endif
3379 int colInA = colIndex[numCol];
3380 int nbCol = A->blocksize1[colInA];
3381 if(colInA)
3382 nbCol -= A->blocksize1[colInA - 1];
3383 if(numCol)
3384 C->blocksize1[numCol] = C->blocksize1[numCol - 1] + nbCol;
3385 else
3386 C->blocksize1[numCol] = nbCol;
3387 }
3388 #ifdef SBM_DEBUG_SBM_COL_PERM
3389 titi = fopen("SBM_column_permutation_output.txt", "w");
3390 SBM_write_in_fileForScilab(C, titi);
3391 fclose(titi);
3392 #endif
3393 }
3394