Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
blas_bsr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 
15 /*---------------------------------*/
16 /*-- Public Functions --*/
17 /*---------------------------------*/
18 
31  const REAL alpha)
32 {
33  const INT nnz = A->NNZ;
34  const INT nb = A->nb;
35 
36  // A direct calculation can be written as:
37  fasp_blas_array_ax(nnz*nb*nb, alpha, A->val);
38 }
39 
59 void fasp_blas_dbsr_aAxpby (const REAL alpha,
60  dBSRmat *A,
61  REAL *x,
62  const REAL beta,
63  REAL *y )
64 {
65  /* members of A */
66  INT ROW = A->ROW;
67  INT nb = A->nb;
68  INT *IA = A->IA;
69  INT *JA = A->JA;
70  REAL *val = A->val;
71 
72  /* local variables */
73  INT size = ROW*nb;
74  INT jump = nb*nb;
75  INT i,j,k,iend;
76  REAL temp;
77  REAL *pA = NULL;
78  REAL *px0 = NULL;
79  REAL *py0 = NULL;
80  REAL *py = NULL;
81 
82  INT nthreads = 1, use_openmp = FALSE;
83 
84 #ifdef _OPENMP
85  if ( ROW > OPENMP_HOLDS ) {
86  use_openmp = TRUE;
87  nthreads = FASP_GET_NUM_THREADS();
88  }
89 #endif
90 
91  //----------------------------------------------
92  // Treat (alpha == 0.0) computation
93  //----------------------------------------------
94 
95  if (alpha == 0.0) {
96  fasp_blas_array_ax(size, beta, y);
97  return;
98  }
99 
100  //-------------------------------------------------
101  // y = (beta/alpha)*y
102  //-------------------------------------------------
103 
104  temp = beta / alpha;
105  if (temp != 1.0) {
106  if (temp == 0.0) {
107  memset(y, 0X0, size*sizeof(REAL));
108  }
109  else {
110  //for (i = size; i--; ) y[i] *= temp; // modified by Xiaozhe, 03/11/2011
111  fasp_blas_array_ax(size, temp, y);
112  }
113  }
114 
115  //-----------------------------------------------------------------
116  // y += A*x (Core Computation)
117  // each non-zero block elements are stored in row-major order
118  //-----------------------------------------------------------------
119 
120  switch (nb)
121  {
122  case 2:
123  {
124  if (use_openmp) {
125  INT myid, mybegin, myend;
126 #ifdef _OPENMP
127 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend)
128 #endif
129  for (myid =0; myid < nthreads; myid++) {
130  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
131  for (i=mybegin; i < myend; ++i) {
132  py0 = &y[i*2];
133  iend = IA[i+1];
134  for (k = IA[i]; k < iend; ++k) {
135  j = JA[k];
136  pA = val+k*4; // &val[k*jump];
137  px0 = x+j*2; // &x[j*nb];
138  py = py0;
139  fasp_blas_smat_ypAx_nc2( pA, px0, py );
140  }
141  }
142  }
143  }
144  else {
145  for (i = 0; i < ROW; ++i) {
146  py0 = &y[i*2];
147  iend = IA[i+1];
148  for (k = IA[i]; k < iend; ++k) {
149  j = JA[k];
150  pA = val+k*4; // &val[k*jump];
151  px0 = x+j*2; // &x[j*nb];
152  py = py0;
153  fasp_blas_smat_ypAx_nc2( pA, px0, py );
154  }
155  }
156  }
157  }
158  break;
159  case 3:
160  {
161  if (use_openmp) {
162  INT myid, mybegin, myend;
163 #ifdef _OPENMP
164 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
165 #endif
166  for (myid =0; myid < nthreads; myid++) {
167  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
168  for (i=mybegin; i < myend; ++i) {
169  py0 = &y[i*3];
170  iend = IA[i+1];
171  for (k = IA[i]; k < iend; ++k) {
172  j = JA[k];
173  pA = val+k*9; // &val[k*jump];
174  px0 = x+j*3; // &x[j*nb];
175  py = py0;
176  fasp_blas_smat_ypAx_nc3( pA, px0, py );
177  }
178  }
179  }
180  }
181  else {
182  for (i = 0; i < ROW; ++i) {
183  py0 = &y[i*3];
184  iend = IA[i+1];
185  for (k = IA[i]; k < iend; ++k) {
186  j = JA[k];
187  pA = val+k*9; // &val[k*jump];
188  px0 = x+j*3; // &x[j*nb];
189  py = py0;
190  fasp_blas_smat_ypAx_nc3( pA, px0, py );
191  }
192  }
193  }
194  }
195  break;
196 
197  case 5:
198  {
199  if (use_openmp) {
200  INT myid, mybegin, myend;
201 #ifdef _OPENMP
202 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
203 #endif
204  for (myid =0; myid < nthreads; myid++) {
205  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
206  for (i=mybegin; i < myend; ++i) {
207  py0 = &y[i*5];
208  iend = IA[i+1];
209  for (k = IA[i]; k < iend; ++k) {
210  j = JA[k];
211  pA = val+k*25; // &val[k*jump];
212  px0 = x+j*5; // &x[j*nb];
213  py = py0;
214  fasp_blas_smat_ypAx_nc5( pA, px0, py );
215  }
216  }
217  }
218  }
219  else {
220  for (i = 0; i < ROW; ++i) {
221  py0 = &y[i*5];
222  iend = IA[i+1];
223  for (k = IA[i]; k < iend; ++k) {
224  j = JA[k];
225  pA = val+k*25; // &val[k*jump];
226  px0 = x+j*5; // &x[j*nb];
227  py = py0;
228  fasp_blas_smat_ypAx_nc5( pA, px0, py );
229  }
230  }
231  }
232  }
233  break;
234  case 7:
235  {
236  if (use_openmp) {
237  INT myid, mybegin, myend;
238 #ifdef _OPENMP
239 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
240 #endif
241  for (myid =0; myid < nthreads; myid++) {
242  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
243  for (i=mybegin; i < myend; ++i) {
244  py0 = &y[i*7];
245  iend = IA[i+1];
246  for (k = IA[i]; k < iend; ++k) {
247  j = JA[k];
248  pA = val+k*49; // &val[k*jump];
249  px0 = x+j*7; // &x[j*nb];
250  py = py0;
251  fasp_blas_smat_ypAx_nc7( pA, px0, py );
252  }
253  }
254  }
255  }
256  else {
257  for (i = 0; i < ROW; ++i) {
258  py0 = &y[i*7];
259  iend = IA[i+1];
260  for (k = IA[i]; k < iend; ++k) {
261  j = JA[k];
262  pA = val+k*49; // &val[k*jump];
263  px0 = x+j*7; // &x[j*nb];
264  py = py0;
265  fasp_blas_smat_ypAx_nc7( pA, px0, py );
266  }
267  }
268  }
269  }
270  break;
271 
272  default:
273  {
274  if (use_openmp) {
275  INT myid, mybegin, myend;
276 #ifdef _OPENMP
277 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend)
278 #endif
279  for (myid =0; myid < nthreads; myid++) {
280  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
281  for (i=mybegin; i < myend; ++i) {
282  py0 = &y[i*nb];
283  iend = IA[i+1];
284  for (k = IA[i]; k < iend; ++k) {
285  j = JA[k];
286  pA = val+k*jump; // &val[k*jump];
287  px0 = x+j*nb; // &x[j*nb];
288  py = py0;
289  fasp_blas_smat_ypAx( pA, px0, py, nb );
290  }
291  }
292  }
293  }
294  else {
295  for (i = 0; i < ROW; ++i) {
296  py0 = &y[i*nb];
297  iend = IA[i+1];
298  for (k = IA[i]; k < iend; ++k) {
299  j = JA[k];
300  pA = val+k*jump; // &val[k*jump];
301  px0 = x+j*nb; // &x[j*nb];
302  py = py0;
303  fasp_blas_smat_ypAx( pA, px0, py, nb );
304  }
305  }
306  }
307  }
308  break;
309  }
310 
311  //------------------------------------------
312  // y = alpha*y
313  //------------------------------------------
314 
315  if (alpha != 1.0) {
316  fasp_blas_array_ax(size, alpha, y);
317  }
318 }
319 
337 void fasp_blas_dbsr_aAxpy (const REAL alpha,
338  dBSRmat *A,
339  REAL *x,
340  REAL *y)
341 {
342  /* members of A */
343  INT ROW = A->ROW;
344  INT nb = A->nb;
345  INT *IA = A->IA;
346  INT *JA = A->JA;
347  REAL *val = A->val;
348 
349  /* local variables */
350  INT size = ROW*nb;
351  INT jump = nb*nb;
352  INT i,j,k, iend;
353  REAL temp = 0.0;
354  REAL *pA = NULL;
355  REAL *px0 = NULL;
356  REAL *py0 = NULL;
357  REAL *py = NULL;
358 
359  INT nthreads = 1, use_openmp = FALSE;
360 
361 #ifdef _OPENMP
362  if ( ROW > OPENMP_HOLDS ) {
363  use_openmp = TRUE;
364  nthreads = FASP_GET_NUM_THREADS();
365  }
366 #endif
367 
368  //----------------------------------------------
369  // Treat (alpha == 0.0) computation
370  //----------------------------------------------
371 
372  if (alpha == 0.0){
373  return; // Nothing to compute
374  }
375 
376  //-------------------------------------------------
377  // y = (1.0/alpha)*y
378  //-------------------------------------------------
379 
380  if (alpha != 1.0){
381  temp = 1.0 / alpha;
382  fasp_blas_array_ax(size, temp, y);
383  }
384 
385  //-----------------------------------------------------------------
386  // y += A*x (Core Computation)
387  // each non-zero block elements are stored in row-major order
388  //-----------------------------------------------------------------
389 
390  switch (nb)
391  {
392  case 2:
393  {
394  if (use_openmp) {
395  INT myid, mybegin, myend;
396 #ifdef _OPENMP
397 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend)
398 #endif
399  for (myid =0; myid < nthreads; myid++) {
400  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
401  for (i=mybegin; i < myend; ++i) {
402  py0 = &y[i*2];
403  iend = IA[i+1];
404  for (k = IA[i]; k < iend; ++k) {
405  j = JA[k];
406  pA = val+k*4; // &val[k*jump];
407  px0 = x+j*2; // &x[j*nb];
408  py = py0;
409  fasp_blas_smat_ypAx_nc2( pA, px0, py );
410  }
411  }
412  }
413  }
414  else {
415  for (i = 0; i < ROW; ++i) {
416  py0 = &y[i*2];
417  iend = IA[i+1];
418  for (k = IA[i]; k < iend; ++k) {
419  j = JA[k];
420  pA = val+k*4; // &val[k*jump];
421  px0 = x+j*2; // &x[j*nb];
422  py = py0;
423  fasp_blas_smat_ypAx_nc2( pA, px0, py );
424  }
425  }
426  }
427  }
428  break;
429 
430  case 3:
431  {
432  if (use_openmp) {
433  INT myid, mybegin, myend;
434 #ifdef _OPENMP
435 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
436 #endif
437  for (myid =0; myid < nthreads; myid++) {
438  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
439  for (i=mybegin; i < myend; ++i) {
440  py0 = &y[i*3];
441  iend = IA[i+1];
442  for (k = IA[i]; k < iend; ++k) {
443  j = JA[k];
444  pA = val+k*9; // &val[k*jump];
445  px0 = x+j*3; // &x[j*nb];
446  py = py0;
447  fasp_blas_smat_ypAx_nc3( pA, px0, py );
448  }
449  }
450  }
451  }
452  else {
453  for (i = 0; i < ROW; ++i){
454  py0 = &y[i*3];
455  iend = IA[i+1];
456  for (k = IA[i]; k < iend; ++k) {
457  j = JA[k];
458  pA = val+k*9; // &val[k*jump];
459  px0 = x+j*3; // &x[j*nb];
460  py = py0;
461  fasp_blas_smat_ypAx_nc3( pA, px0, py );
462  }
463  }
464  }
465  }
466  break;
467 
468  case 5:
469  {
470  if (use_openmp) {
471  INT myid, mybegin, myend;
472 #ifdef _OPENMP
473 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
474 #endif
475  for (myid =0; myid < nthreads; myid++) {
476  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
477  for (i=mybegin; i < myend; ++i) {
478  py0 = &y[i*5];
479  iend = IA[i+1];
480  for (k = IA[i]; k < iend; ++k) {
481  j = JA[k];
482  pA = val+k*25; // &val[k*jump];
483  px0 = x+j*5; // &x[j*nb];
484  py = py0;
485  fasp_blas_smat_ypAx_nc5( pA, px0, py );
486  }
487  }
488  }
489  }
490  else {
491  for (i = 0; i < ROW; ++i){
492  py0 = &y[i*5];
493  iend = IA[i+1];
494  for (k = IA[i]; k < iend; ++k) {
495  j = JA[k];
496  pA = val+k*25; // &val[k*jump];
497  px0 = x+j*5; // &x[j*nb];
498  py = py0;
499  fasp_blas_smat_ypAx_nc5( pA, px0, py );
500  }
501  }
502  }
503  }
504  break;
505  case 7:
506  {
507  if (use_openmp) {
508  INT myid, mybegin, myend;
509 #ifdef _OPENMP
510 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend )
511 #endif
512  for (myid =0; myid < nthreads; myid++) {
513  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
514  for (i=mybegin; i < myend; ++i) {
515  py0 = &y[i*7];
516  iend = IA[i+1];
517  for (k = IA[i]; k < iend; ++k) {
518  j = JA[k];
519  pA = val+k*49; // &val[k*jump];
520  px0 = x+j*7; // &x[j*nb];
521  py = py0;
522  fasp_blas_smat_ypAx_nc7( pA, px0, py );
523  }
524  }
525  }
526  }
527  else {
528  for (i = 0; i < ROW; ++i) {
529  py0 = &y[i*7];
530  iend = IA[i+1];
531  for (k = IA[i]; k < iend; ++k) {
532  j = JA[k];
533  pA = val+k*49; // &val[k*jump];
534  px0 = x+j*7; // &x[j*nb];
535  py = py0;
536  fasp_blas_smat_ypAx_nc7( pA, px0, py );
537  }
538 
539  }
540  }
541  }
542  break;
543 
544  default:
545  {
546  if (use_openmp) {
547  INT myid, mybegin, myend;
548 #ifdef _OPENMP
549 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, pA, px0, py,iend)
550 #endif
551  for (myid =0; myid < nthreads; myid++) {
552  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
553  for (i=mybegin; i < myend; ++i) {
554  py0 = &y[i*nb];
555  iend = IA[i+1];
556  for (k = IA[i]; k < iend; ++k) {
557  j = JA[k];
558  pA = val+k*jump; // &val[k*jump];
559  px0 = x+j*nb; // &x[j*nb];
560  py = py0;
561  fasp_blas_smat_ypAx( pA, px0, py, nb );
562  }
563 
564  }
565  }
566  }
567  else {
568  for (i = 0; i < ROW; ++i) {
569  py0 = &y[i*nb];
570  iend = IA[i+1];
571  for (k = IA[i]; k < iend; ++k) {
572  j = JA[k];
573  pA = val+k*jump; // &val[k*jump];
574  px0 = x+j*nb; // &x[j*nb];
575  py = py0;
576  fasp_blas_smat_ypAx( pA, px0, py, nb );
577  }
578 
579  }
580  }
581  }
582  break;
583  }
584 
585  //------------------------------------------
586  // y = alpha*y
587  //------------------------------------------
588 
589  if (alpha != 1.0){
590  fasp_blas_array_ax(size, alpha, y);
591  }
592  return;
593 }
594 
610 void fasp_blas_dbsr_aAxpy_agg (const REAL alpha,
611  dBSRmat *A,
612  REAL *x,
613  REAL *y)
614 {
615  /* members of A */
616  INT ROW = A->ROW;
617  INT nb = A->nb;
618  INT *IA = A->IA;
619  INT *JA = A->JA;
620 
621  /* local variables */
622  INT size = ROW*nb;
623  INT i,j,k, iend;
624  REAL temp = 0.0;
625  REAL *px0 = NULL, *py0 = NULL, *py = NULL;
626 
627  INT nthreads = 1, use_openmp = FALSE;
628 
629 #ifdef _OPENMP
630  if ( ROW > OPENMP_HOLDS ) {
631  use_openmp = TRUE;
632  nthreads = FASP_GET_NUM_THREADS();
633  }
634 #endif
635 
636  //----------------------------------------------
637  // Treat (alpha == 0.0) computation
638  //----------------------------------------------
639 
640  if (alpha == 0.0){
641  return; // Nothing to compute
642  }
643 
644  //-------------------------------------------------
645  // y = (1.0/alpha)*y
646  //-------------------------------------------------
647 
648  if (alpha != 1.0){
649  temp = 1.0 / alpha;
650  fasp_blas_array_ax(size, temp, y);
651  }
652 
653  //-----------------------------------------------------------------
654  // y += A*x (Core Computation)
655  // each non-zero block elements are stored in row-major order
656  //-----------------------------------------------------------------
657 
658  switch (nb)
659  {
660  case 2:
661  {
662  if (use_openmp) {
663  INT myid, mybegin, myend;
664 #ifdef _OPENMP
665 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, px0, py,iend)
666 #endif
667  for (myid =0; myid < nthreads; myid++) {
668  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
669  for (i=mybegin; i < myend; ++i) {
670  py0 = &y[i*2];
671  iend = IA[i+1];
672  for (k = IA[i]; k < iend; ++k) {
673  j = JA[k];
674  px0 = x+j*2; // &x[j*nb];
675  py = py0;
676  py[0] += px0[0];
677  py[1] += px0[1];
678  }
679  }
680  }
681  }
682  else {
683  for (i = 0; i < ROW; ++i) {
684  py0 = &y[i*2];
685  iend = IA[i+1];
686  for (k = IA[i]; k < iend; ++k) {
687  j = JA[k];
688  px0 = x+j*2; // &x[j*nb];
689  py = py0;
690  py[0] += px0[0];
691  py[1] += px0[1];
692  }
693  }
694  }
695  }
696  break;
697 
698  case 3:
699  {
700  if (use_openmp) {
701  INT myid, mybegin, myend;
702 #ifdef _OPENMP
703 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, px0, py,iend )
704 #endif
705  for (myid =0; myid < nthreads; myid++) {
706  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
707  for (i=mybegin; i < myend; ++i) {
708  py0 = &y[i*3];
709  iend = IA[i+1];
710  for (k = IA[i]; k < iend; ++k) {
711  j = JA[k];
712  px0 = x+j*3; // &x[j*nb];
713  py = py0;
714  py[0] += px0[0];
715  py[1] += px0[1];
716  py[2] += px0[2];
717  }
718  }
719  }
720  }
721  else {
722  for (i = 0; i < ROW; ++i){
723  py0 = &y[i*3];
724  iend = IA[i+1];
725  for (k = IA[i]; k < iend; ++k) {
726  j = JA[k];
727  px0 = x+j*3; // &x[j*nb];
728  py = py0;
729  py[0] += px0[0];
730  py[1] += px0[1];
731  py[2] += px0[2];
732  }
733  }
734  }
735  }
736  break;
737 
738  case 5:
739  {
740  if (use_openmp) {
741  INT myid, mybegin, myend;
742 #ifdef _OPENMP
743 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, px0, py,iend )
744 #endif
745  for (myid =0; myid < nthreads; myid++) {
746  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
747  for (i=mybegin; i < myend; ++i) {
748  py0 = &y[i*5];
749  iend = IA[i+1];
750  for (k = IA[i]; k < iend; ++k) {
751  j = JA[k];
752  px0 = x+j*5; // &x[j*nb];
753  py = py0;
754  py[0] += px0[0];
755  py[1] += px0[1];
756  py[2] += px0[2];
757  py[3] += px0[3];
758  py[4] += px0[4];
759  }
760  }
761  }
762  }
763  else {
764  for (i = 0; i < ROW; ++i){
765  py0 = &y[i*5];
766  iend = IA[i+1];
767  for (k = IA[i]; k < iend; ++k) {
768  j = JA[k];
769  px0 = x+j*5; // &x[j*nb];
770  py = py0;
771  py[0] += px0[0];
772  py[1] += px0[1];
773  py[2] += px0[2];
774  py[3] += px0[3];
775  py[4] += px0[4];
776  }
777  }
778  }
779  }
780  break;
781 
782  case 7:
783  {
784  if (use_openmp) {
785  INT myid, mybegin, myend;
786 #ifdef _OPENMP
787 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, px0, py,iend )
788 #endif
789  for (myid =0; myid < nthreads; myid++) {
790  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
791  for (i=mybegin; i < myend; ++i) {
792  py0 = &y[i*7];
793  iend = IA[i+1];
794  for (k = IA[i]; k < iend; ++k) {
795  j = JA[k];
796  px0 = x+j*7; // &x[j*nb];
797  py = py0;
798  py[0] += px0[0];
799  py[1] += px0[1];
800  py[2] += px0[2];
801  py[3] += px0[3];
802  py[4] += px0[4];
803  py[5] += px0[5];
804  py[6] += px0[6];
805  }
806  }
807  }
808  }
809  else {
810  for (i = 0; i < ROW; ++i) {
811  py0 = &y[i*7];
812  iend = IA[i+1];
813  for (k = IA[i]; k < iend; ++k) {
814  j = JA[k];
815  px0 = x+j*7; // &x[j*nb];
816  py = py0;
817  py[0] += px0[0];
818  py[1] += px0[1];
819  py[2] += px0[2];
820  py[3] += px0[3];
821  py[4] += px0[4];
822  py[5] += px0[5];
823  py[6] += px0[6];
824  }
825 
826  }
827  }
828  }
829  break;
830 
831  default:
832  {
833  if (use_openmp) {
834  INT myid, mybegin, myend;
835 #ifdef _OPENMP
836 #pragma omp parallel for private(myid, mybegin, myend, i, py0, k, j, px0, py,iend)
837 #endif
838  for (myid =0; myid < nthreads; myid++) {
839  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
840  for (i=mybegin; i < myend; ++i) {
841  py0 = &y[i*nb];
842  iend = IA[i+1];
843  for (k = IA[i]; k < iend; ++k) {
844  j = JA[k];
845  px0 = x+j*nb; // &x[j*nb];
846  py = py0;
847  fasp_blas_array_axpy(nb, 1.0, px0, py);
848  }
849 
850  }
851  }
852  }
853  else {
854  for (i = 0; i < ROW; ++i) {
855  py0 = &y[i*nb];
856  iend = IA[i+1];
857  for (k = IA[i]; k < iend; ++k) {
858  j = JA[k];
859  px0 = x+j*nb; // &x[j*nb];
860  py = py0;
861  fasp_blas_array_axpy(nb, 1.0, px0, py);
862  }
863 
864  }
865  }
866  }
867  break;
868  }
869 
870  //------------------------------------------
871  // y = alpha*y
872  //------------------------------------------
873 
874  if ( alpha != 1.0 ) fasp_blas_array_ax(size, alpha, y);
875 
876  return;
877 }
878 
896  REAL *x,
897  REAL *y)
898 {
899  /* members of A */
900  INT ROW = A->ROW;
901  INT nb = A->nb;
902  INT *IA = A->IA;
903  INT *JA = A->JA;
904  REAL *val = A->val;
905 
906  /* local variables */
907  INT size = ROW*nb;
908  INT jump = nb*nb;
909  INT i,j,k, num_nnz_row;
910  REAL *pA = NULL;
911  REAL *px0 = NULL;
912  REAL *py0 = NULL;
913  REAL *py = NULL;
914 
915  INT use_openmp = FALSE;
916 
917 #ifdef _OPENMP
918  INT myid, mybegin, myend, nthreads;
919  if ( ROW > OPENMP_HOLDS ) {
920  use_openmp = TRUE;
921  nthreads = FASP_GET_NUM_THREADS();
922  }
923 #endif
924 
925  //-----------------------------------------------------------------
926  // zero out 'y'
927  //-----------------------------------------------------------------
928  fasp_array_set(size, y, 0.0);
929 
930  //-----------------------------------------------------------------
931  // y = A*x (Core Computation)
932  // each non-zero block elements are stored in row-major order
933  //-----------------------------------------------------------------
934 
935  switch (nb)
936  {
937  case 3:
938  {
939  if (use_openmp) {
940 #ifdef _OPENMP
941 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
942  {
943  myid = omp_get_thread_num();
944  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
945  for (i=mybegin; i < myend; ++i)
946  {
947  py0 = &y[i*3];
948  num_nnz_row = IA[i+1] - IA[i];
949  switch(num_nnz_row)
950  {
951  case 3:
952  k = IA[i];
953  j = JA[k];
954  pA = val+k*9;
955  px0 = x+j*3;
956  py = py0;
957  fasp_blas_smat_ypAx_nc3( pA, px0, py );
958 
959  k ++;
960  j = JA[k];
961  pA = val+k*9;
962  px0 = x+j*3;
963  py = py0;
964  fasp_blas_smat_ypAx_nc3( pA, px0, py );
965 
966  k ++;
967  j = JA[k];
968  pA = val+k*9;
969  px0 = x+j*3;
970  py = py0;
971  fasp_blas_smat_ypAx_nc3( pA, px0, py );
972 
973  break;
974  case 4:
975  k = IA[i];
976  j = JA[k];
977  pA = val+k*9;
978  px0 = x+j*3;
979  py = py0;
980  fasp_blas_smat_ypAx_nc3( pA, px0, py );
981 
982  k ++;
983  j = JA[k];
984  pA = val+k*9;
985  px0 = x+j*3;
986  py = py0;
987  fasp_blas_smat_ypAx_nc3( pA, px0, py );
988 
989  k ++;
990  j = JA[k];
991  pA = val+k*9;
992  px0 = x+j*3;
993  py = py0;
994  fasp_blas_smat_ypAx_nc3( pA, px0, py );
995 
996  k ++;
997  j = JA[k];
998  pA = val+k*9;
999  px0 = x+j*3;
1000  py = py0;
1001  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1002 
1003  break;
1004  case 5:
1005  k = IA[i];
1006  j = JA[k];
1007  pA = val+k*9;
1008  px0 = x+j*3;
1009  py = py0;
1010  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1011 
1012  k ++;
1013  j = JA[k];
1014  pA = val+k*9;
1015  px0 = x+j*3;
1016  py = py0;
1017  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1018 
1019  k ++;
1020  j = JA[k];
1021  pA = val+k*9;
1022  px0 = x+j*3;
1023  py = py0;
1024  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1025 
1026  k ++;
1027  j = JA[k];
1028  pA = val+k*9;
1029  px0 = x+j*3;
1030  py = py0;
1031  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1032 
1033  k ++;
1034  j = JA[k];
1035  pA = val+k*9;
1036  px0 = x+j*3;
1037  py = py0;
1038  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1039 
1040  break;
1041  case 6:
1042  k = IA[i];
1043  j = JA[k];
1044  pA = val+k*9;
1045  px0 = x+j*3;
1046  py = py0;
1047  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1048 
1049  k ++;
1050  j = JA[k];
1051  pA = val+k*9;
1052  px0 = x+j*3;
1053  py = py0;
1054  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1055 
1056  k ++;
1057  j = JA[k];
1058  pA = val+k*9;
1059  px0 = x+j*3;
1060  py = py0;
1061  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1062 
1063  k ++;
1064  j = JA[k];
1065  pA = val+k*9;
1066  px0 = x+j*3;
1067  py = py0;
1068  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1069 
1070  k ++;
1071  j = JA[k];
1072  pA = val+k*9;
1073  px0 = x+j*3;
1074  py = py0;
1075  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1076 
1077  k ++;
1078  j = JA[k];
1079  pA = val+k*9;
1080  px0 = x+j*3;
1081  py = py0;
1082  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1083 
1084  break;
1085  case 7:
1086  k = IA[i];
1087  j = JA[k];
1088  pA = val+k*9;
1089  px0 = x+j*3;
1090  py = py0;
1091  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1092 
1093  k ++;
1094  j = JA[k];
1095  pA = val+k*9;
1096  px0 = x+j*3;
1097  py = py0;
1098  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1099 
1100  k ++;
1101  j = JA[k];
1102  pA = val+k*9;
1103  px0 = x+j*3;
1104  py = py0;
1105  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1106 
1107  k ++;
1108  j = JA[k];
1109  pA = val+k*9;
1110  px0 = x+j*3;
1111  py = py0;
1112  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1113 
1114  k ++;
1115  j = JA[k];
1116  pA = val+k*9;
1117  px0 = x+j*3;
1118  py = py0;
1119  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1120 
1121  k ++;
1122  j = JA[k];
1123  pA = val+k*9;
1124  px0 = x+j*3;
1125  py = py0;
1126  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1127 
1128  k ++;
1129  j = JA[k];
1130  pA = val+k*9;
1131  px0 = x+j*3;
1132  py = py0;
1133  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1134 
1135  break;
1136  default:
1137  for (k = IA[i]; k < IA[i+1]; ++k)
1138  {
1139  j = JA[k];
1140  pA = val+k*9;
1141  px0 = x+j*3;
1142  py = py0;
1143  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1144  }
1145  break;
1146  }
1147  }
1148  }
1149 #endif
1150  }
1151  else {
1152  for (i = 0; i < ROW; ++i)
1153  {
1154  py0 = &y[i*3];
1155  num_nnz_row = IA[i+1] - IA[i];
1156  switch(num_nnz_row)
1157  {
1158  case 3:
1159  k = IA[i];
1160  j = JA[k];
1161  pA = val+k*9; // &val[k*jump];
1162  px0 = x+j*3; // &x[j*nb];
1163  py = py0;
1164  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1165 
1166  k ++;
1167  j = JA[k];
1168  pA = val+k*9; // &val[k*jump];
1169  px0 = x+j*3; // &x[j*nb];
1170  py = py0;
1171  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1172 
1173  k ++;
1174  j = JA[k];
1175  pA = val+k*9; // &val[k*jump];
1176  px0 = x+j*3; // &x[j*nb];
1177  py = py0;
1178  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1179 
1180  break;
1181  case 4:
1182  k = IA[i];
1183  j = JA[k];
1184  pA = val+k*9; // &val[k*jump];
1185  px0 = x+j*3; // &x[j*nb];
1186  py = py0;
1187  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1188 
1189  k ++;
1190  j = JA[k];
1191  pA = val+k*9; // &val[k*jump];
1192  px0 = x+j*3; // &x[j*nb];
1193  py = py0;
1194  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1195 
1196  k ++;
1197  j = JA[k];
1198  pA = val+k*9; // &val[k*jump];
1199  px0 = x+j*3; // &x[j*nb];
1200  py = py0;
1201  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1202 
1203  k ++;
1204  j = JA[k];
1205  pA = val+k*9; // &val[k*jump];
1206  px0 = x+j*3; // &x[j*nb];
1207  py = py0;
1208  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1209 
1210  break;
1211  case 5:
1212  k = IA[i];
1213  j = JA[k];
1214  pA = val+k*9; // &val[k*jump];
1215  px0 = x+j*3; // &x[j*nb];
1216  py = py0;
1217  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1218 
1219  k ++;
1220  j = JA[k];
1221  pA = val+k*9; // &val[k*jump];
1222  px0 = x+j*3; // &x[j*nb];
1223  py = py0;
1224  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1225 
1226  k ++;
1227  j = JA[k];
1228  pA = val+k*9; // &val[k*jump];
1229  px0 = x+j*3; // &x[j*nb];
1230  py = py0;
1231  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1232 
1233  k ++;
1234  j = JA[k];
1235  pA = val+k*9; // &val[k*jump];
1236  px0 = x+j*3; // &x[j*nb];
1237  py = py0;
1238  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1239 
1240  k ++;
1241  j = JA[k];
1242  pA = val+k*9; // &val[k*jump];
1243  px0 = x+j*3; // &x[j*nb];
1244  py = py0;
1245  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1246 
1247  break;
1248  case 6:
1249  k = IA[i];
1250  j = JA[k];
1251  pA = val+k*9; // &val[k*jump];
1252  px0 = x+j*3; // &x[j*nb];
1253  py = py0;
1254  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1255 
1256  k ++;
1257  j = JA[k];
1258  pA = val+k*9; // &val[k*jump];
1259  px0 = x+j*3; // &x[j*nb];
1260  py = py0;
1261  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1262 
1263  k ++;
1264  j = JA[k];
1265  pA = val+k*9; // &val[k*jump];
1266  px0 = x+j*3; // &x[j*nb];
1267  py = py0;
1268  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1269 
1270  k ++;
1271  j = JA[k];
1272  pA = val+k*9; // &val[k*jump];
1273  px0 = x+j*3; // &x[j*nb];
1274  py = py0;
1275  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1276 
1277  k ++;
1278  j = JA[k];
1279  pA = val+k*9; // &val[k*jump];
1280  px0 = x+j*3; // &x[j*nb];
1281  py = py0;
1282  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1283 
1284  k ++;
1285  j = JA[k];
1286  pA = val+k*9; // &val[k*jump];
1287  px0 = x+j*3; // &x[j*nb];
1288  py = py0;
1289  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1290 
1291  break;
1292  case 7:
1293  k = IA[i];
1294  j = JA[k];
1295  pA = val+k*9; // &val[k*jump];
1296  px0 = x+j*3; // &x[j*nb];
1297  py = py0;
1298  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1299 
1300  k ++;
1301  j = JA[k];
1302  pA = val+k*9; // &val[k*jump];
1303  px0 = x+j*3; // &x[j*nb];
1304  py = py0;
1305  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1306 
1307  k ++;
1308  j = JA[k];
1309  pA = val+k*9; // &val[k*jump];
1310  px0 = x+j*3; // &x[j*nb];
1311  py = py0;
1312  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1313 
1314  k ++;
1315  j = JA[k];
1316  pA = val+k*9; // &val[k*jump];
1317  px0 = x+j*3; // &x[j*nb];
1318  py = py0;
1319  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1320 
1321  k ++;
1322  j = JA[k];
1323  pA = val+k*9; // &val[k*jump];
1324  px0 = x+j*3; // &x[j*nb];
1325  py = py0;
1326  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1327 
1328  k ++;
1329  j = JA[k];
1330  pA = val+k*9; // &val[k*jump];
1331  px0 = x+j*3; // &x[j*nb];
1332  py = py0;
1333  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1334 
1335  k ++;
1336  j = JA[k];
1337  pA = val+k*9; // &val[k*jump];
1338  px0 = x+j*3; // &x[j*nb];
1339  py = py0;
1340  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1341 
1342  break;
1343  default:
1344  for (k = IA[i]; k < IA[i+1]; ++k)
1345  {
1346  j = JA[k];
1347  pA = val+k*9; // &val[k*jump];
1348  px0 = x+j*3; // &x[j*nb];
1349  py = py0;
1350  fasp_blas_smat_ypAx_nc3( pA, px0, py );
1351  }
1352  break;
1353  }
1354  }
1355  }
1356  }
1357  break;
1358 
1359  case 5:
1360  {
1361  if (use_openmp) {
1362 #ifdef _OPENMP
1363 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
1364  {
1365  myid = omp_get_thread_num();
1366  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1367  for (i=mybegin; i < myend; ++i)
1368  {
1369  py0 = &y[i*5];
1370  num_nnz_row = IA[i+1] - IA[i];
1371  switch(num_nnz_row)
1372  {
1373  case 3:
1374  k = IA[i];
1375  j = JA[k];
1376  pA = val+k*25; // &val[k*jump];
1377  px0 = x+j*5; // &x[j*nb];
1378  py = py0;
1379  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1380 
1381  k ++;
1382  j = JA[k];
1383  pA = val+k*25; // &val[k*jump];
1384  px0 = x+j*5; // &x[j*nb];
1385  py = py0;
1386  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1387 
1388  k ++;
1389  j = JA[k];
1390  pA = val+k*25; // &val[k*jump];
1391  px0 = x+j*5; // &x[j*nb];
1392  py = py0;
1393  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1394 
1395  break;
1396  case 4:
1397  k = IA[i];
1398  j = JA[k];
1399  pA = val+k*25; // &val[k*jump];
1400  px0 = x+j*5; // &x[j*nb];
1401  py = py0;
1402  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1403 
1404  k ++;
1405  j = JA[k];
1406  pA = val+k*25; // &val[k*jump];
1407  px0 = x+j*5; // &x[j*nb];
1408  py = py0;
1409  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1410 
1411  k ++;
1412  j = JA[k];
1413  pA = val+k*25; // &val[k*jump];
1414  px0 = x+j*5; // &x[j*nb];
1415  py = py0;
1416  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1417 
1418  k ++;
1419  j = JA[k];
1420  pA = val+k*25; // &val[k*jump];
1421  px0 = x+j*5; // &x[j*nb];
1422  py = py0;
1423  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1424 
1425  break;
1426  case 5:
1427  k = IA[i];
1428  j = JA[k];
1429  pA = val+k*25; // &val[k*jump];
1430  px0 = x+j*5; // &x[j*nb];
1431  py = py0;
1432  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1433 
1434  k ++;
1435  j = JA[k];
1436  pA = val+k*25; // &val[k*jump];
1437  px0 = x+j*5; // &x[j*nb];
1438  py = py0;
1439  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1440 
1441  k ++;
1442  j = JA[k];
1443  pA = val+k*25; // &val[k*jump];
1444  px0 = x+j*5; // &x[j*nb];
1445  py = py0;
1446  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1447 
1448  k ++;
1449  j = JA[k];
1450  pA = val+k*25; // &val[k*jump];
1451  px0 = x+j*5; // &x[j*nb];
1452  py = py0;
1453  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1454 
1455  k ++;
1456  j = JA[k];
1457  pA = val+k*25; // &val[k*jump];
1458  px0 = x+j*5; // &x[j*nb];
1459  py = py0;
1460  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1461 
1462  break;
1463  case 6:
1464  k = IA[i];
1465  j = JA[k];
1466  pA = val+k*25; // &val[k*jump];
1467  px0 = x+j*5; // &x[j*nb];
1468  py = py0;
1469  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1470 
1471  k ++;
1472  j = JA[k];
1473  pA = val+k*25; // &val[k*jump];
1474  px0 = x+j*5; // &x[j*nb];
1475  py = py0;
1476  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1477 
1478  k ++;
1479  j = JA[k];
1480  pA = val+k*25; // &val[k*jump];
1481  px0 = x+j*5; // &x[j*nb];
1482  py = py0;
1483  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1484 
1485  k ++;
1486  j = JA[k];
1487  pA = val+k*25; // &val[k*jump];
1488  px0 = x+j*5; // &x[j*nb];
1489  py = py0;
1490  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1491 
1492  k ++;
1493  j = JA[k];
1494  pA = val+k*25; // &val[k*jump];
1495  px0 = x+j*5; // &x[j*nb];
1496  py = py0;
1497  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1498 
1499  k ++;
1500  j = JA[k];
1501  pA = val+k*25; // &val[k*jump];
1502  px0 = x+j*5; // &x[j*nb];
1503  py = py0;
1504  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1505 
1506  break;
1507  case 7:
1508  k = IA[i];
1509  j = JA[k];
1510  pA = val+k*25; // &val[k*jump];
1511  px0 = x+j*5; // &x[j*nb];
1512  py = py0;
1513  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1514 
1515  k ++;
1516  j = JA[k];
1517  pA = val+k*25; // &val[k*jump];
1518  px0 = x+j*5; // &x[j*nb];
1519  py = py0;
1520  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1521 
1522  k ++;
1523  j = JA[k];
1524  pA = val+k*25; // &val[k*jump];
1525  px0 = x+j*5; // &x[j*nb];
1526  py = py0;
1527  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1528 
1529  k ++;
1530  j = JA[k];
1531  pA = val+k*25; // &val[k*jump];
1532  px0 = x+j*5; // &x[j*nb];
1533  py = py0;
1534  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1535 
1536  k ++;
1537  j = JA[k];
1538  pA = val+k*25; // &val[k*jump];
1539  px0 = x+j*5; // &x[j*nb];
1540  py = py0;
1541  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1542 
1543  k ++;
1544  j = JA[k];
1545  pA = val+k*25; // &val[k*jump];
1546  px0 = x+j*5; // &x[j*nb];
1547  py = py0;
1548  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1549 
1550  k ++;
1551  j = JA[k];
1552  pA = val+k*25; // &val[k*jump];
1553  px0 = x+j*5; // &x[j*nb];
1554  py = py0;
1555  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1556 
1557  break;
1558  default:
1559  for (k = IA[i]; k < IA[i+1]; ++k)
1560  {
1561  j = JA[k];
1562  pA = val+k*25; // &val[k*jump];
1563  px0 = x+j*5; // &x[j*nb];
1564  py = py0;
1565  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1566  }
1567  break;
1568  }
1569  }
1570  }
1571 #endif
1572  }
1573  else {
1574  for (i = 0; i < ROW; ++i)
1575  {
1576  py0 = &y[i*5];
1577  num_nnz_row = IA[i+1] - IA[i];
1578  switch(num_nnz_row)
1579  {
1580  case 3:
1581  k = IA[i];
1582  j = JA[k];
1583  pA = val+k*25; // &val[k*jump];
1584  px0 = x+j*5; // &x[j*nb];
1585  py = py0;
1586  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1587 
1588  k ++;
1589  j = JA[k];
1590  pA = val+k*25; // &val[k*jump];
1591  px0 = x+j*5; // &x[j*nb];
1592  py = py0;
1593  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1594 
1595  k ++;
1596  j = JA[k];
1597  pA = val+k*25; // &val[k*jump];
1598  px0 = x+j*5; // &x[j*nb];
1599  py = py0;
1600  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1601 
1602  break;
1603  case 4:
1604  k = IA[i];
1605  j = JA[k];
1606  pA = val+k*25; // &val[k*jump];
1607  px0 = x+j*5; // &x[j*nb];
1608  py = py0;
1609  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1610 
1611  k ++;
1612  j = JA[k];
1613  pA = val+k*25; // &val[k*jump];
1614  px0 = x+j*5; // &x[j*nb];
1615  py = py0;
1616  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1617 
1618  k ++;
1619  j = JA[k];
1620  pA = val+k*25; // &val[k*jump];
1621  px0 = x+j*5; // &x[j*nb];
1622  py = py0;
1623  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1624 
1625  k ++;
1626  j = JA[k];
1627  pA = val+k*25; // &val[k*jump];
1628  px0 = x+j*5; // &x[j*nb];
1629  py = py0;
1630  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1631 
1632  break;
1633  case 5:
1634  k = IA[i];
1635  j = JA[k];
1636  pA = val+k*25; // &val[k*jump];
1637  px0 = x+j*5; // &x[j*nb];
1638  py = py0;
1639  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1640 
1641  k ++;
1642  j = JA[k];
1643  pA = val+k*25; // &val[k*jump];
1644  px0 = x+j*5; // &x[j*nb];
1645  py = py0;
1646  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1647 
1648  k ++;
1649  j = JA[k];
1650  pA = val+k*25; // &val[k*jump];
1651  px0 = x+j*5; // &x[j*nb];
1652  py = py0;
1653  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1654 
1655  k ++;
1656  j = JA[k];
1657  pA = val+k*25; // &val[k*jump];
1658  px0 = x+j*5; // &x[j*nb];
1659  py = py0;
1660  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1661 
1662  k ++;
1663  j = JA[k];
1664  pA = val+k*25; // &val[k*jump];
1665  px0 = x+j*5; // &x[j*nb];
1666  py = py0;
1667  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1668 
1669  break;
1670  case 6:
1671  k = IA[i];
1672  j = JA[k];
1673  pA = val+k*25; // &val[k*jump];
1674  px0 = x+j*5; // &x[j*nb];
1675  py = py0;
1676  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1677 
1678  k ++;
1679  j = JA[k];
1680  pA = val+k*25; // &val[k*jump];
1681  px0 = x+j*5; // &x[j*nb];
1682  py = py0;
1683  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1684 
1685  k ++;
1686  j = JA[k];
1687  pA = val+k*25; // &val[k*jump];
1688  px0 = x+j*5; // &x[j*nb];
1689  py = py0;
1690  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1691 
1692  k ++;
1693  j = JA[k];
1694  pA = val+k*25; // &val[k*jump];
1695  px0 = x+j*5; // &x[j*nb];
1696  py = py0;
1697  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1698 
1699  k ++;
1700  j = JA[k];
1701  pA = val+k*25; // &val[k*jump];
1702  px0 = x+j*5; // &x[j*nb];
1703  py = py0;
1704  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1705 
1706  k ++;
1707  j = JA[k];
1708  pA = val+k*25; // &val[k*jump];
1709  px0 = x+j*5; // &x[j*nb];
1710  py = py0;
1711  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1712 
1713  break;
1714  case 7:
1715  k = IA[i];
1716  j = JA[k];
1717  pA = val+k*25; // &val[k*jump];
1718  px0 = x+j*5; // &x[j*nb];
1719  py = py0;
1720  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1721 
1722  k ++;
1723  j = JA[k];
1724  pA = val+k*25; // &val[k*jump];
1725  px0 = x+j*5; // &x[j*nb];
1726  py = py0;
1727  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1728 
1729  k ++;
1730  j = JA[k];
1731  pA = val+k*25; // &val[k*jump];
1732  px0 = x+j*5; // &x[j*nb];
1733  py = py0;
1734  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1735 
1736  k ++;
1737  j = JA[k];
1738  pA = val+k*25; // &val[k*jump];
1739  px0 = x+j*5; // &x[j*nb];
1740  py = py0;
1741  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1742 
1743  k ++;
1744  j = JA[k];
1745  pA = val+k*25; // &val[k*jump];
1746  px0 = x+j*5; // &x[j*nb];
1747  py = py0;
1748  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1749 
1750  k ++;
1751  j = JA[k];
1752  pA = val+k*25; // &val[k*jump];
1753  px0 = x+j*5; // &x[j*nb];
1754  py = py0;
1755  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1756 
1757  k ++;
1758  j = JA[k];
1759  pA = val+k*25; // &val[k*jump];
1760  px0 = x+j*5; // &x[j*nb];
1761  py = py0;
1762  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1763 
1764  break;
1765  default:
1766  for (k = IA[i]; k < IA[i+1]; ++k)
1767  {
1768  j = JA[k];
1769  pA = val+k*25; // &val[k*jump];
1770  px0 = x+j*5; // &x[j*nb];
1771  py = py0;
1772  fasp_blas_smat_ypAx_nc5( pA, px0, py );
1773  }
1774  break;
1775  }
1776  }
1777  }
1778  }
1779  break;
1780 
1781  case 7:
1782  {
1783  if (use_openmp) {
1784 #ifdef _OPENMP
1785 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
1786  {
1787  myid = omp_get_thread_num();
1788  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
1789  for (i=mybegin; i < myend; ++i)
1790  {
1791  py0 = &y[i*7];
1792  num_nnz_row = IA[i+1] - IA[i];
1793  switch(num_nnz_row)
1794  {
1795  case 3:
1796  k = IA[i];
1797  j = JA[k];
1798  pA = val+k*49; // &val[k*jump];
1799  px0 = x+j*7; // &x[j*nb];
1800  py = py0;
1801  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1802 
1803  k ++;
1804  j = JA[k];
1805  pA = val+k*49; // &val[k*jump];
1806  px0 = x+j*7; // &x[j*nb];
1807  py = py0;
1808  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1809 
1810  k ++;
1811  j = JA[k];
1812  pA = val+k*49; // &val[k*jump];
1813  px0 = x+j*7; // &x[j*nb];
1814  py = py0;
1815  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1816 
1817  break;
1818  case 4:
1819  k = IA[i];
1820  j = JA[k];
1821  pA = val+k*49; // &val[k*jump];
1822  px0 = x+j*7; // &x[j*nb];
1823  py = py0;
1824  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1825 
1826  k ++;
1827  j = JA[k];
1828  pA = val+k*49; // &val[k*jump];
1829  px0 = x+j*7; // &x[j*nb];
1830  py = py0;
1831  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1832 
1833  k ++;
1834  j = JA[k];
1835  pA = val+k*49; // &val[k*jump];
1836  px0 = x+j*7; // &x[j*nb];
1837  py = py0;
1838  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1839 
1840  k ++;
1841  j = JA[k];
1842  pA = val+k*49; // &val[k*jump];
1843  px0 = x+j*7; // &x[j*nb];
1844  py = py0;
1845  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1846 
1847  break;
1848  case 5:
1849  k = IA[i];
1850  j = JA[k];
1851  pA = val+k*49; // &val[k*jump];
1852  px0 = x+j*7; // &x[j*nb];
1853  py = py0;
1854  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1855 
1856  k ++;
1857  j = JA[k];
1858  pA = val+k*49; // &val[k*jump];
1859  px0 = x+j*7; // &x[j*nb];
1860  py = py0;
1861  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1862 
1863  k ++;
1864  j = JA[k];
1865  pA = val+k*49; // &val[k*jump];
1866  px0 = x+j*7; // &x[j*nb];
1867  py = py0;
1868  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1869 
1870  k ++;
1871  j = JA[k];
1872  pA = val+k*49; // &val[k*jump];
1873  px0 = x+j*7; // &x[j*nb];
1874  py = py0;
1875  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1876 
1877  k ++;
1878  j = JA[k];
1879  pA = val+k*49; // &val[k*jump];
1880  px0 = x+j*7; // &x[j*nb];
1881  py = py0;
1882  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1883 
1884  break;
1885  case 6:
1886  k = IA[i];
1887  j = JA[k];
1888  pA = val+k*49; // &val[k*jump];
1889  px0 = x+j*7; // &x[j*nb];
1890  py = py0;
1891  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1892 
1893  k ++;
1894  j = JA[k];
1895  pA = val+k*49; // &val[k*jump];
1896  px0 = x+j*7; // &x[j*nb];
1897  py = py0;
1898  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1899 
1900  k ++;
1901  j = JA[k];
1902  pA = val+k*49; // &val[k*jump];
1903  px0 = x+j*7; // &x[j*nb];
1904  py = py0;
1905  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1906 
1907  k ++;
1908  j = JA[k];
1909  pA = val+k*49; // &val[k*jump];
1910  px0 = x+j*7; // &x[j*nb];
1911  py = py0;
1912  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1913 
1914  k ++;
1915  j = JA[k];
1916  pA = val+k*49; // &val[k*jump];
1917  px0 = x+j*7; // &x[j*nb];
1918  py = py0;
1919  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1920 
1921  k ++;
1922  j = JA[k];
1923  pA = val+k*49; // &val[k*jump];
1924  px0 = x+j*7; // &x[j*nb];
1925  py = py0;
1926  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1927 
1928  break;
1929  case 7:
1930  k = IA[i];
1931  j = JA[k];
1932  pA = val+k*49; // &val[k*jump];
1933  px0 = x+j*7; // &x[j*nb];
1934  py = py0;
1935  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1936 
1937  k ++;
1938  j = JA[k];
1939  pA = val+k*49; // &val[k*jump];
1940  px0 = x+j*7; // &x[j*nb];
1941  py = py0;
1942  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1943 
1944  k ++;
1945  j = JA[k];
1946  pA = val+k*49; // &val[k*jump];
1947  px0 = x+j*7; // &x[j*nb];
1948  py = py0;
1949  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1950 
1951  k ++;
1952  j = JA[k];
1953  pA = val+k*49; // &val[k*jump];
1954  px0 = x+j*7; // &x[j*nb];
1955  py = py0;
1956  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1957 
1958  k ++;
1959  j = JA[k];
1960  pA = val+k*49; // &val[k*jump];
1961  px0 = x+j*7; // &x[j*nb];
1962  py = py0;
1963  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1964 
1965  k ++;
1966  j = JA[k];
1967  pA = val+k*49; // &val[k*jump];
1968  px0 = x+j*7; // &x[j*nb];
1969  py = py0;
1970  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1971 
1972  k ++;
1973  j = JA[k];
1974  pA = val+k*49; // &val[k*jump];
1975  px0 = x+j*7; // &x[j*nb];
1976  py = py0;
1977  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1978 
1979  break;
1980  default:
1981  for (k = IA[i]; k < IA[i+1]; ++k)
1982  {
1983  j = JA[k];
1984  pA = val+k*49; // &val[k*jump];
1985  px0 = x+j*7; // &x[j*nb];
1986  py = py0;
1987  fasp_blas_smat_ypAx_nc7( pA, px0, py );
1988  }
1989  break;
1990  }
1991  }
1992  }
1993 #endif
1994  }
1995  else {
1996  for (i = 0; i < ROW; ++i)
1997  {
1998  py0 = &y[i*7];
1999  num_nnz_row = IA[i+1] - IA[i];
2000  switch(num_nnz_row)
2001  {
2002  case 3:
2003  k = IA[i];
2004  j = JA[k];
2005  pA = val+k*49; // &val[k*jump];
2006  px0 = x+j*7; // &x[j*nb];
2007  py = py0;
2008  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2009 
2010  k ++;
2011  j = JA[k];
2012  pA = val+k*49; // &val[k*jump];
2013  px0 = x+j*7; // &x[j*nb];
2014  py = py0;
2015  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2016 
2017  k ++;
2018  j = JA[k];
2019  pA = val+k*49; // &val[k*jump];
2020  px0 = x+j*7; // &x[j*nb];
2021  py = py0;
2022  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2023 
2024  break;
2025  case 4:
2026  k = IA[i];
2027  j = JA[k];
2028  pA = val+k*49; // &val[k*jump];
2029  px0 = x+j*7; // &x[j*nb];
2030  py = py0;
2031  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2032 
2033  k ++;
2034  j = JA[k];
2035  pA = val+k*49; // &val[k*jump];
2036  px0 = x+j*7; // &x[j*nb];
2037  py = py0;
2038  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2039 
2040  k ++;
2041  j = JA[k];
2042  pA = val+k*49; // &val[k*jump];
2043  px0 = x+j*7; // &x[j*nb];
2044  py = py0;
2045  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2046 
2047  k ++;
2048  j = JA[k];
2049  pA = val+k*49; // &val[k*jump];
2050  px0 = x+j*7; // &x[j*nb];
2051  py = py0;
2052  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2053 
2054  break;
2055  case 5:
2056  k = IA[i];
2057  j = JA[k];
2058  pA = val+k*49; // &val[k*jump];
2059  px0 = x+j*7; // &x[j*nb];
2060  py = py0;
2061  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2062 
2063  k ++;
2064  j = JA[k];
2065  pA = val+k*49; // &val[k*jump];
2066  px0 = x+j*7; // &x[j*nb];
2067  py = py0;
2068  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2069 
2070  k ++;
2071  j = JA[k];
2072  pA = val+k*49; // &val[k*jump];
2073  px0 = x+j*7; // &x[j*nb];
2074  py = py0;
2075  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2076 
2077  k ++;
2078  j = JA[k];
2079  pA = val+k*49; // &val[k*jump];
2080  px0 = x+j*7; // &x[j*nb];
2081  py = py0;
2082  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2083 
2084  k ++;
2085  j = JA[k];
2086  pA = val+k*49; // &val[k*jump];
2087  px0 = x+j*7; // &x[j*nb];
2088  py = py0;
2089  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2090 
2091  break;
2092  case 6:
2093  k = IA[i];
2094  j = JA[k];
2095  pA = val+k*49; // &val[k*jump];
2096  px0 = x+j*7; // &x[j*nb];
2097  py = py0;
2098  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2099 
2100  k ++;
2101  j = JA[k];
2102  pA = val+k*49; // &val[k*jump];
2103  px0 = x+j*7; // &x[j*nb];
2104  py = py0;
2105  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2106 
2107  k ++;
2108  j = JA[k];
2109  pA = val+k*49; // &val[k*jump];
2110  px0 = x+j*7; // &x[j*nb];
2111  py = py0;
2112  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2113 
2114  k ++;
2115  j = JA[k];
2116  pA = val+k*49; // &val[k*jump];
2117  px0 = x+j*7; // &x[j*nb];
2118  py = py0;
2119  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2120 
2121  k ++;
2122  j = JA[k];
2123  pA = val+k*49; // &val[k*jump];
2124  px0 = x+j*7; // &x[j*nb];
2125  py = py0;
2126  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2127 
2128  k ++;
2129  j = JA[k];
2130  pA = val+k*49; // &val[k*jump];
2131  px0 = x+j*7; // &x[j*nb];
2132  py = py0;
2133  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2134 
2135  break;
2136  case 7:
2137  k = IA[i];
2138  j = JA[k];
2139  pA = val+k*49; // &val[k*jump];
2140  px0 = x+j*7; // &x[j*nb];
2141  py = py0;
2142  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2143 
2144  k ++;
2145  j = JA[k];
2146  pA = val+k*49; // &val[k*jump];
2147  px0 = x+j*7; // &x[j*nb];
2148  py = py0;
2149  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2150 
2151  k ++;
2152  j = JA[k];
2153  pA = val+k*49; // &val[k*jump];
2154  px0 = x+j*7; // &x[j*nb];
2155  py = py0;
2156  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2157 
2158  k ++;
2159  j = JA[k];
2160  pA = val+k*49; // &val[k*jump];
2161  px0 = x+j*7; // &x[j*nb];
2162  py = py0;
2163  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2164 
2165  k ++;
2166  j = JA[k];
2167  pA = val+k*49; // &val[k*jump];
2168  px0 = x+j*7; // &x[j*nb];
2169  py = py0;
2170  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2171 
2172  k ++;
2173  j = JA[k];
2174  pA = val+k*49; // &val[k*jump];
2175  px0 = x+j*7; // &x[j*nb];
2176  py = py0;
2177  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2178 
2179  k ++;
2180  j = JA[k];
2181  pA = val+k*49; // &val[k*jump];
2182  px0 = x+j*7; // &x[j*nb];
2183  py = py0;
2184  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2185 
2186  break;
2187  default:
2188  for (k = IA[i]; k < IA[i+1]; ++k)
2189  {
2190  j = JA[k];
2191  pA = val+k*49; // &val[k*jump];
2192  px0 = x+j*7; // &x[j*nb];
2193  py = py0;
2194  fasp_blas_smat_ypAx_nc7( pA, px0, py );
2195  }
2196  break;
2197  }
2198  }
2199  }
2200  }
2201  break;
2202 
2203  default:
2204  {
2205  if (use_openmp) {
2206 #ifdef _OPENMP
2207 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
2208  {
2209  myid = omp_get_thread_num();
2210  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
2211  for (i=mybegin; i < myend; ++i)
2212  {
2213  py0 = &y[i*nb];
2214  num_nnz_row = IA[i+1] - IA[i];
2215  switch(num_nnz_row)
2216  {
2217  case 3:
2218  k = IA[i];
2219  j = JA[k];
2220  pA = val+k*jump; // &val[k*jump];
2221  px0 = x+j*nb; // &x[j*nb];
2222  py = py0;
2223  fasp_blas_smat_ypAx( pA, px0, py, nb );
2224 
2225  k ++;
2226  j = JA[k];
2227  pA = val+k*jump; // &val[k*jump];
2228  px0 = x+j*nb; // &x[j*nb];
2229  py = py0;
2230  fasp_blas_smat_ypAx( pA, px0, py, nb );
2231 
2232  k ++;
2233  j = JA[k];
2234  pA = val+k*jump; // &val[k*jump];
2235  px0 = x+j*nb; // &x[j*nb];
2236  py = py0;
2237  fasp_blas_smat_ypAx( pA, px0, py, nb );
2238 
2239  break;
2240  case 4:
2241  k = IA[i];
2242  j = JA[k];
2243  pA = val+k*jump; // &val[k*jump];
2244  px0 = x+j*nb; // &x[j*nb];
2245  py = py0;
2246  fasp_blas_smat_ypAx( pA, px0, py, nb );
2247 
2248  k ++;
2249  j = JA[k];
2250  pA = val+k*jump; // &val[k*jump];
2251  px0 = x+j*nb; // &x[j*nb];
2252  py = py0;
2253  fasp_blas_smat_ypAx( pA, px0, py, nb );
2254 
2255  k ++;
2256  j = JA[k];
2257  pA = val+k*jump; // &val[k*jump];
2258  px0 = x+j*nb; // &x[j*nb];
2259  py = py0;
2260  fasp_blas_smat_ypAx( pA, px0, py, nb );
2261 
2262  k ++;
2263  j = JA[k];
2264  pA = val+k*jump; // &val[k*jump];
2265  px0 = x+j*nb; // &x[j*nb];
2266  py = py0;
2267  fasp_blas_smat_ypAx( pA, px0, py, nb );
2268 
2269  break;
2270  case 5:
2271  k = IA[i];
2272  j = JA[k];
2273  pA = val+k*jump; // &val[k*jump];
2274  px0 = x+j*nb; // &x[j*nb];
2275  py = py0;
2276  fasp_blas_smat_ypAx( pA, px0, py, nb );
2277 
2278  k ++;
2279  j = JA[k];
2280  pA = val+k*jump; // &val[k*jump];
2281  px0 = x+j*nb; // &x[j*nb];
2282  py = py0;
2283  fasp_blas_smat_ypAx( pA, px0, py, nb );
2284 
2285  k ++;
2286  j = JA[k];
2287  pA = val+k*jump; // &val[k*jump];
2288  px0 = x+j*nb; // &x[j*nb];
2289  py = py0;
2290  fasp_blas_smat_ypAx( pA, px0, py, nb );
2291 
2292  k ++;
2293  j = JA[k];
2294  pA = val+k*jump; // &val[k*jump];
2295  px0 = x+j*nb; // &x[j*nb];
2296  py = py0;
2297  fasp_blas_smat_ypAx( pA, px0, py, nb );
2298 
2299  k ++;
2300  j = JA[k];
2301  pA = val+k*jump; // &val[k*jump];
2302  px0 = x+j*nb; // &x[j*nb];
2303  py = py0;
2304  fasp_blas_smat_ypAx( pA, px0, py, nb );
2305 
2306  break;
2307  case 6:
2308  k = IA[i];
2309  j = JA[k];
2310  pA = val+k*jump; // &val[k*jump];
2311  px0 = x+j*nb; // &x[j*nb];
2312  py = py0;
2313  fasp_blas_smat_ypAx( pA, px0, py, nb );
2314 
2315  k ++;
2316  j = JA[k];
2317  pA = val+k*jump; // &val[k*jump];
2318  px0 = x+j*nb; // &x[j*nb];
2319  py = py0;
2320  fasp_blas_smat_ypAx( pA, px0, py, nb );
2321 
2322  k ++;
2323  j = JA[k];
2324  pA = val+k*jump; // &val[k*jump];
2325  px0 = x+j*nb; // &x[j*nb];
2326  py = py0;
2327  fasp_blas_smat_ypAx( pA, px0, py, nb );
2328 
2329  k ++;
2330  j = JA[k];
2331  pA = val+k*jump; // &val[k*jump];
2332  px0 = x+j*nb; // &x[j*nb];
2333  py = py0;
2334  fasp_blas_smat_ypAx( pA, px0, py, nb );
2335 
2336  k ++;
2337  j = JA[k];
2338  pA = val+k*jump; // &val[k*jump];
2339  px0 = x+j*nb; // &x[j*nb];
2340  py = py0;
2341  fasp_blas_smat_ypAx( pA, px0, py, nb );
2342 
2343  k ++;
2344  j = JA[k];
2345  pA = val+k*jump; // &val[k*jump];
2346  px0 = x+j*nb; // &x[j*nb];
2347  py = py0;
2348  fasp_blas_smat_ypAx( pA, px0, py, nb );
2349 
2350  break;
2351  case 7:
2352  k = IA[i];
2353  j = JA[k];
2354  pA = val+k*jump; // &val[k*jump];
2355  px0 = x+j*nb; // &x[j*nb];
2356  py = py0;
2357  fasp_blas_smat_ypAx( pA, px0, py, nb );
2358 
2359  k ++;
2360  j = JA[k];
2361  pA = val+k*jump; // &val[k*jump];
2362  px0 = x+j*nb; // &x[j*nb];
2363  py = py0;
2364  fasp_blas_smat_ypAx( pA, px0, py, nb );
2365 
2366  k ++;
2367  j = JA[k];
2368  pA = val+k*jump; // &val[k*jump];
2369  px0 = x+j*nb; // &x[j*nb];
2370  py = py0;
2371  fasp_blas_smat_ypAx( pA, px0, py, nb );
2372 
2373  k ++;
2374  j = JA[k];
2375  pA = val+k*jump; // &val[k*jump];
2376  px0 = x+j*nb; // &x[j*nb];
2377  py = py0;
2378  fasp_blas_smat_ypAx( pA, px0, py, nb );
2379 
2380  k ++;
2381  j = JA[k];
2382  pA = val+k*jump; // &val[k*jump];
2383  px0 = x+j*nb; // &x[j*nb];
2384  py = py0;
2385  fasp_blas_smat_ypAx( pA, px0, py, nb );
2386 
2387  k ++;
2388  j = JA[k];
2389  pA = val+k*jump; // &val[k*jump];
2390  px0 = x+j*nb; // &x[j*nb];
2391  py = py0;
2392  fasp_blas_smat_ypAx( pA, px0, py, nb );
2393 
2394  k ++;
2395  j = JA[k];
2396  pA = val+k*jump; // &val[k*jump];
2397  px0 = x+j*nb; // &x[j*nb];
2398  py = py0;
2399  fasp_blas_smat_ypAx( pA, px0, py, nb );
2400 
2401  break;
2402  default:
2403  for (k = IA[i]; k < IA[i+1]; ++k)
2404  {
2405  j = JA[k];
2406  pA = val+k*jump; // &val[k*jump];
2407  px0 = x+j*nb; // &x[j*nb];
2408  py = py0;
2409  fasp_blas_smat_ypAx( pA, px0, py, nb );
2410  }
2411  break;
2412  }
2413  }
2414  }
2415 #endif
2416  }
2417  else {
2418  for (i = 0; i < ROW; ++i)
2419  {
2420  py0 = &y[i*nb];
2421  num_nnz_row = IA[i+1] - IA[i];
2422  switch(num_nnz_row)
2423  {
2424  case 3:
2425  k = IA[i];
2426  j = JA[k];
2427  pA = val+k*jump; // &val[k*jump];
2428  px0 = x+j*nb; // &x[j*nb];
2429  py = py0;
2430  fasp_blas_smat_ypAx( pA, px0, py, nb );
2431 
2432  k ++;
2433  j = JA[k];
2434  pA = val+k*jump; // &val[k*jump];
2435  px0 = x+j*nb; // &x[j*nb];
2436  py = py0;
2437  fasp_blas_smat_ypAx( pA, px0, py, nb );
2438 
2439  k ++;
2440  j = JA[k];
2441  pA = val+k*jump; // &val[k*jump];
2442  px0 = x+j*nb; // &x[j*nb];
2443  py = py0;
2444  fasp_blas_smat_ypAx( pA, px0, py, nb );
2445 
2446  break;
2447  case 4:
2448  k = IA[i];
2449  j = JA[k];
2450  pA = val+k*jump; // &val[k*jump];
2451  px0 = x+j*nb; // &x[j*nb];
2452  py = py0;
2453  fasp_blas_smat_ypAx( pA, px0, py, nb );
2454 
2455  k ++;
2456  j = JA[k];
2457  pA = val+k*jump; // &val[k*jump];
2458  px0 = x+j*nb; // &x[j*nb];
2459  py = py0;
2460  fasp_blas_smat_ypAx( pA, px0, py, nb );
2461 
2462  k ++;
2463  j = JA[k];
2464  pA = val+k*jump; // &val[k*jump];
2465  px0 = x+j*nb; // &x[j*nb];
2466  py = py0;
2467  fasp_blas_smat_ypAx( pA, px0, py, nb );
2468 
2469  k ++;
2470  j = JA[k];
2471  pA = val+k*jump; // &val[k*jump];
2472  px0 = x+j*nb; // &x[j*nb];
2473  py = py0;
2474  fasp_blas_smat_ypAx( pA, px0, py, nb );
2475 
2476  break;
2477  case 5:
2478  k = IA[i];
2479  j = JA[k];
2480  pA = val+k*jump; // &val[k*jump];
2481  px0 = x+j*nb; // &x[j*nb];
2482  py = py0;
2483  fasp_blas_smat_ypAx( pA, px0, py, nb );
2484 
2485  k ++;
2486  j = JA[k];
2487  pA = val+k*jump; // &val[k*jump];
2488  px0 = x+j*nb; // &x[j*nb];
2489  py = py0;
2490  fasp_blas_smat_ypAx( pA, px0, py, nb );
2491 
2492  k ++;
2493  j = JA[k];
2494  pA = val+k*jump; // &val[k*jump];
2495  px0 = x+j*nb; // &x[j*nb];
2496  py = py0;
2497  fasp_blas_smat_ypAx( pA, px0, py, nb );
2498 
2499  k ++;
2500  j = JA[k];
2501  pA = val+k*jump; // &val[k*jump];
2502  px0 = x+j*nb; // &x[j*nb];
2503  py = py0;
2504  fasp_blas_smat_ypAx( pA, px0, py, nb );
2505 
2506  k ++;
2507  j = JA[k];
2508  pA = val+k*jump; // &val[k*jump];
2509  px0 = x+j*nb; // &x[j*nb];
2510  py = py0;
2511  fasp_blas_smat_ypAx( pA, px0, py, nb );
2512 
2513  break;
2514  case 6:
2515  k = IA[i];
2516  j = JA[k];
2517  pA = val+k*jump; // &val[k*jump];
2518  px0 = x+j*nb; // &x[j*nb];
2519  py = py0;
2520  fasp_blas_smat_ypAx( pA, px0, py, nb );
2521 
2522  k ++;
2523  j = JA[k];
2524  pA = val+k*jump; // &val[k*jump];
2525  px0 = x+j*nb; // &x[j*nb];
2526  py = py0;
2527  fasp_blas_smat_ypAx( pA, px0, py, nb );
2528 
2529  k ++;
2530  j = JA[k];
2531  pA = val+k*jump; // &val[k*jump];
2532  px0 = x+j*nb; // &x[j*nb];
2533  py = py0;
2534  fasp_blas_smat_ypAx( pA, px0, py, nb );
2535 
2536  k ++;
2537  j = JA[k];
2538  pA = val+k*jump; // &val[k*jump];
2539  px0 = x+j*nb; // &x[j*nb];
2540  py = py0;
2541  fasp_blas_smat_ypAx( pA, px0, py, nb );
2542 
2543  k ++;
2544  j = JA[k];
2545  pA = val+k*jump; // &val[k*jump];
2546  px0 = x+j*nb; // &x[j*nb];
2547  py = py0;
2548  fasp_blas_smat_ypAx( pA, px0, py, nb );
2549 
2550  k ++;
2551  j = JA[k];
2552  pA = val+k*jump; // &val[k*jump];
2553  px0 = x+j*nb; // &x[j*nb];
2554  py = py0;
2555  fasp_blas_smat_ypAx( pA, px0, py, nb );
2556 
2557  break;
2558  case 7:
2559  k = IA[i];
2560  j = JA[k];
2561  pA = val+k*jump; // &val[k*jump];
2562  px0 = x+j*nb; // &x[j*nb];
2563  py = py0;
2564  fasp_blas_smat_ypAx( pA, px0, py, nb );
2565 
2566  k ++;
2567  j = JA[k];
2568  pA = val+k*jump; // &val[k*jump];
2569  px0 = x+j*nb; // &x[j*nb];
2570  py = py0;
2571  fasp_blas_smat_ypAx( pA, px0, py, nb );
2572 
2573  k ++;
2574  j = JA[k];
2575  pA = val+k*jump; // &val[k*jump];
2576  px0 = x+j*nb; // &x[j*nb];
2577  py = py0;
2578  fasp_blas_smat_ypAx( pA, px0, py, nb );
2579 
2580  k ++;
2581  j = JA[k];
2582  pA = val+k*jump; // &val[k*jump];
2583  px0 = x+j*nb; // &x[j*nb];
2584  py = py0;
2585  fasp_blas_smat_ypAx( pA, px0, py, nb );
2586 
2587  k ++;
2588  j = JA[k];
2589  pA = val+k*jump; // &val[k*jump];
2590  px0 = x+j*nb; // &x[j*nb];
2591  py = py0;
2592  fasp_blas_smat_ypAx( pA, px0, py, nb );
2593 
2594  k ++;
2595  j = JA[k];
2596  pA = val+k*jump; // &val[k*jump];
2597  px0 = x+j*nb; // &x[j*nb];
2598  py = py0;
2599  fasp_blas_smat_ypAx( pA, px0, py, nb );
2600 
2601  k ++;
2602  j = JA[k];
2603  pA = val+k*jump; // &val[k*jump];
2604  px0 = x+j*nb; // &x[j*nb];
2605  py = py0;
2606  fasp_blas_smat_ypAx( pA, px0, py, nb );
2607 
2608  break;
2609  default:
2610  for (k = IA[i]; k < IA[i+1]; ++k)
2611  {
2612  j = JA[k];
2613  pA = val+k*jump; // &val[k*jump];
2614  px0 = x+j*nb; // &x[j*nb];
2615  py = py0;
2616  fasp_blas_smat_ypAx( pA, px0, py, nb );
2617  }
2618  break;
2619  }
2620  }
2621  }
2622  }
2623  break;
2624  }
2625 }
2626 
2642  REAL *x,
2643  REAL *y)
2644 {
2645  /* members of A */
2646  const INT ROW = A->ROW;
2647  const INT nb = A->nb;
2648  const INT size = ROW*nb;
2649 
2650  INT *IA = A->IA;
2651  INT *JA = A->JA;
2652 
2653  /* local variables */
2654  INT i,j,k, num_nnz_row;
2655  REAL *px0 = NULL, *py0 = NULL, *py = NULL;
2656  INT use_openmp = FALSE;
2657 
2658 #ifdef _OPENMP
2659  REAL *val = A->val;
2660  REAL *pA;
2661  INT myid, mybegin, myend, nthreads;
2662  if ( ROW > OPENMP_HOLDS ) {
2663  use_openmp = TRUE;
2664  nthreads = FASP_GET_NUM_THREADS();
2665  }
2666 #endif
2667 
2668  //-----------------------------------------------------------------
2669  // zero out 'y'
2670  //-----------------------------------------------------------------
2671  fasp_array_set(size, y, 0.0);
2672 
2673  //-----------------------------------------------------------------
2674  // y = A*x (Core Computation)
2675  // each non-zero block elements are stored in row-major order
2676  //-----------------------------------------------------------------
2677 
2678  switch (nb)
2679  {
2680  case 3:
2681  {
2682  if (use_openmp) {
2683 #ifdef _OPENMP
2684 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
2685  {
2686  myid = omp_get_thread_num();
2687  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
2688  for (i=mybegin; i < myend; ++i) {
2689  py0 = &y[i*3];
2690  num_nnz_row = IA[i+1] - IA[i];
2691  switch(num_nnz_row) {
2692  case 3:
2693  k = IA[i];
2694  j = JA[k];
2695  pA = val+k*9;
2696  px0 = x+j*3;
2697  py = py0;
2698  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2699 
2700  k ++;
2701  j = JA[k];
2702  pA = val+k*9;
2703  px0 = x+j*3;
2704  py = py0;
2705  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2706 
2707  k ++;
2708  j = JA[k];
2709  pA = val+k*9;
2710  px0 = x+j*3;
2711  py = py0;
2712  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2713 
2714  break;
2715 
2716  case 4:
2717  k = IA[i];
2718  j = JA[k];
2719  pA = val+k*9;
2720  px0 = x+j*3;
2721  py = py0;
2722  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2723 
2724  k ++;
2725  j = JA[k];
2726  pA = val+k*9;
2727  px0 = x+j*3;
2728  py = py0;
2729  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2730 
2731  k ++;
2732  j = JA[k];
2733  pA = val+k*9;
2734  px0 = x+j*3;
2735  py = py0;
2736  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2737 
2738  k ++;
2739  j = JA[k];
2740  pA = val+k*9;
2741  px0 = x+j*3;
2742  py = py0;
2743  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2744 
2745  break;
2746 
2747  case 5:
2748  k = IA[i];
2749  j = JA[k];
2750  pA = val+k*9;
2751  px0 = x+j*3;
2752  py = py0;
2753  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2754 
2755  k ++;
2756  j = JA[k];
2757  pA = val+k*9;
2758  px0 = x+j*3;
2759  py = py0;
2760  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2761 
2762  k ++;
2763  j = JA[k];
2764  pA = val+k*9;
2765  px0 = x+j*3;
2766  py = py0;
2767  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2768 
2769  k ++;
2770  j = JA[k];
2771  pA = val+k*9;
2772  px0 = x+j*3;
2773  py = py0;
2774  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2775 
2776  k ++;
2777  j = JA[k];
2778  pA = val+k*9;
2779  px0 = x+j*3;
2780  py = py0;
2781  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2782 
2783  break;
2784 
2785  case 6:
2786  k = IA[i];
2787  j = JA[k];
2788  pA = val+k*9;
2789  px0 = x+j*3;
2790  py = py0;
2791  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2792 
2793  k ++;
2794  j = JA[k];
2795  pA = val+k*9;
2796  px0 = x+j*3;
2797  py = py0;
2798  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2799 
2800  k ++;
2801  j = JA[k];
2802  pA = val+k*9;
2803  px0 = x+j*3;
2804  py = py0;
2805  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2806 
2807  k ++;
2808  j = JA[k];
2809  pA = val+k*9;
2810  px0 = x+j*3;
2811  py = py0;
2812  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2813 
2814  k ++;
2815  j = JA[k];
2816  pA = val+k*9;
2817  px0 = x+j*3;
2818  py = py0;
2819  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2820 
2821  k ++;
2822  j = JA[k];
2823  pA = val+k*9;
2824  px0 = x+j*3;
2825  py = py0;
2826  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2827 
2828  break;
2829 
2830  case 7:
2831  k = IA[i];
2832  j = JA[k];
2833  pA = val+k*9;
2834  px0 = x+j*3;
2835  py = py0;
2836  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2837 
2838  k ++;
2839  j = JA[k];
2840  pA = val+k*9;
2841  px0 = x+j*3;
2842  py = py0;
2843  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2844 
2845  k ++;
2846  j = JA[k];
2847  pA = val+k*9;
2848  px0 = x+j*3;
2849  py = py0;
2850  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2851 
2852  k ++;
2853  j = JA[k];
2854  pA = val+k*9;
2855  px0 = x+j*3;
2856  py = py0;
2857  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2858 
2859  k ++;
2860  j = JA[k];
2861  pA = val+k*9;
2862  px0 = x+j*3;
2863  py = py0;
2864  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2865 
2866  k ++;
2867  j = JA[k];
2868  pA = val+k*9;
2869  px0 = x+j*3;
2870  py = py0;
2871  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2872 
2873  k ++;
2874  j = JA[k];
2875  pA = val+k*9;
2876  px0 = x+j*3;
2877  py = py0;
2878  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2879 
2880  break;
2881 
2882  default:
2883  for (k = IA[i]; k < IA[i+1]; ++k)
2884  {
2885  j = JA[k];
2886  pA = val+k*9;
2887  px0 = x+j*3;
2888  py = py0;
2889  fasp_blas_smat_ypAx_nc3( pA, px0, py );
2890  }
2891  break;
2892  }
2893  }
2894  }
2895 #endif
2896  }
2897  else {
2898  for (i = 0; i < ROW; ++i) {
2899  py0 = &y[i*3];
2900  num_nnz_row = IA[i+1] - IA[i];
2901  switch(num_nnz_row) {
2902  case 3:
2903  k = IA[i];
2904  j = JA[k];
2905  px0 = x+j*3; // &x[j*nb];
2906  py = py0;
2907  py[0] += px0[0];
2908  py[1] += px0[1];
2909  py[2] += px0[2];
2910 
2911  k ++;
2912  j = JA[k];
2913  px0 = x+j*3; // &x[j*nb];
2914  py = py0;
2915  py[0] += px0[0];
2916  py[1] += px0[1];
2917  py[2] += px0[2];
2918 
2919  k ++;
2920  j = JA[k];
2921  px0 = x+j*3; // &x[j*nb];
2922  py = py0;
2923  py[0] += px0[0];
2924  py[1] += px0[1];
2925  py[2] += px0[2];
2926 
2927  break;
2928 
2929  case 4:
2930  k = IA[i];
2931  j = JA[k];
2932  px0 = x+j*3; // &x[j*nb];
2933  py = py0;
2934  py[0] += px0[0];
2935  py[1] += px0[1];
2936  py[2] += px0[2];
2937 
2938  k ++;
2939  j = JA[k];
2940  px0 = x+j*3; // &x[j*nb];
2941  py = py0;
2942  py[0] += px0[0];
2943  py[1] += px0[1];
2944  py[2] += px0[2];
2945 
2946  k ++;
2947  j = JA[k];
2948  px0 = x+j*3; // &x[j*nb];
2949  py = py0;
2950  py[0] += px0[0];
2951  py[1] += px0[1];
2952  py[2] += px0[2];
2953 
2954  k ++;
2955  j = JA[k];
2956  px0 = x+j*3; // &x[j*nb];
2957  py = py0;
2958  py[0] += px0[0];
2959  py[1] += px0[1];
2960  py[2] += px0[2];
2961 
2962  break;
2963 
2964  case 5:
2965  k = IA[i];
2966  j = JA[k];
2967  px0 = x+j*3; // &x[j*nb];
2968  py = py0;
2969  py[0] += px0[0];
2970  py[1] += px0[1];
2971  py[2] += px0[2];
2972 
2973  k ++;
2974  j = JA[k];
2975  px0 = x+j*3; // &x[j*nb];
2976  py = py0;
2977  py[0] += px0[0];
2978  py[1] += px0[1];
2979  py[2] += px0[2];
2980 
2981  k ++;
2982  j = JA[k];
2983  px0 = x+j*3; // &x[j*nb];
2984  py = py0;
2985  py[0] += px0[0];
2986  py[1] += px0[1];
2987  py[2] += px0[2];
2988 
2989  k ++;
2990  j = JA[k];
2991  px0 = x+j*3; // &x[j*nb];
2992  py = py0;
2993  py[0] += px0[0];
2994  py[1] += px0[1];
2995  py[2] += px0[2];
2996 
2997  k ++;
2998  j = JA[k];
2999  px0 = x+j*3; // &x[j*nb];
3000  py = py0;
3001  py[0] += px0[0];
3002  py[1] += px0[1];
3003  py[2] += px0[2];
3004 
3005  break;
3006 
3007  case 6:
3008  k = IA[i];
3009  j = JA[k];
3010  px0 = x+j*3; // &x[j*nb];
3011  py = py0;
3012  py[0] += px0[0];
3013  py[1] += px0[1];
3014  py[2] += px0[2];
3015 
3016  k ++;
3017  j = JA[k];
3018  px0 = x+j*3; // &x[j*nb];
3019  py = py0;
3020  py[0] += px0[0];
3021  py[1] += px0[1];
3022  py[2] += px0[2];
3023 
3024  k ++;
3025  j = JA[k];
3026  px0 = x+j*3; // &x[j*nb];
3027  py = py0;
3028  py[0] += px0[0];
3029  py[1] += px0[1];
3030  py[2] += px0[2];
3031 
3032  k ++;
3033  j = JA[k];
3034  px0 = x+j*3; // &x[j*nb];
3035  py = py0;
3036  py[0] += px0[0];
3037  py[1] += px0[1];
3038  py[2] += px0[2];
3039 
3040  k ++;
3041  j = JA[k];
3042  px0 = x+j*3; // &x[j*nb];
3043  py = py0;
3044  py[0] += px0[0];
3045  py[1] += px0[1];
3046  py[2] += px0[2];
3047 
3048  k ++;
3049  j = JA[k];
3050  px0 = x+j*3; // &x[j*nb];
3051  py = py0;
3052  py[0] += px0[0];
3053  py[1] += px0[1];
3054  py[2] += px0[2];
3055 
3056  break;
3057 
3058  case 7:
3059  k = IA[i];
3060  j = JA[k];
3061  px0 = x+j*3; // &x[j*nb];
3062  py = py0;
3063  py[0] += px0[0];
3064  py[1] += px0[1];
3065  py[2] += px0[2];
3066 
3067  k ++;
3068  j = JA[k];
3069  px0 = x+j*3; // &x[j*nb];
3070  py = py0;
3071  py[0] += px0[0];
3072  py[1] += px0[1];
3073  py[2] += px0[2];
3074 
3075  k ++;
3076  j = JA[k];
3077  px0 = x+j*3; // &x[j*nb];
3078  py = py0;
3079  py[0] += px0[0];
3080  py[1] += px0[1];
3081  py[2] += px0[2];
3082 
3083  k ++;
3084  j = JA[k];
3085  px0 = x+j*3; // &x[j*nb];
3086  py = py0;
3087  py[0] += px0[0];
3088  py[1] += px0[1];
3089  py[2] += px0[2];
3090 
3091  k ++;
3092  j = JA[k];
3093  px0 = x+j*3; // &x[j*nb];
3094  py = py0;
3095  py[0] += px0[0];
3096  py[1] += px0[1];
3097  py[2] += px0[2];
3098 
3099  k ++;
3100  j = JA[k];
3101  px0 = x+j*3; // &x[j*nb];
3102  py = py0;
3103  py[0] += px0[0];
3104  py[1] += px0[1];
3105  py[2] += px0[2];
3106 
3107  k ++;
3108  j = JA[k];
3109  px0 = x+j*3; // &x[j*nb];
3110  py = py0;
3111  py[0] += px0[0];
3112  py[1] += px0[1];
3113  py[2] += px0[2];
3114 
3115  break;
3116 
3117  default:
3118  for (k = IA[i]; k < IA[i+1]; ++k) {
3119  j = JA[k];
3120  px0 = x+j*3; // &x[j*nb];
3121  py = py0;
3122  py[0] += px0[0];
3123  py[1] += px0[1];
3124  py[2] += px0[2];
3125  }
3126  break;
3127  }
3128  }
3129  }
3130  }
3131  break;
3132 
3133  case 5:
3134  {
3135  if (use_openmp) {
3136 #ifdef _OPENMP
3137 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
3138  {
3139  myid = omp_get_thread_num();
3140  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
3141  for (i=mybegin; i < myend; ++i) {
3142  py0 = &y[i*5];
3143  num_nnz_row = IA[i+1] - IA[i];
3144  switch(num_nnz_row) {
3145 
3146  case 3:
3147  k = IA[i];
3148  j = JA[k];
3149  pA = val+k*25; // &val[k*jump];
3150  px0 = x+j*5; // &x[j*nb];
3151  py = py0;
3152  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3153 
3154  k ++;
3155  j = JA[k];
3156  pA = val+k*25; // &val[k*jump];
3157  px0 = x+j*5; // &x[j*nb];
3158  py = py0;
3159  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3160 
3161  k ++;
3162  j = JA[k];
3163  pA = val+k*25; // &val[k*jump];
3164  px0 = x+j*5; // &x[j*nb];
3165  py = py0;
3166  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3167 
3168  break;
3169 
3170  case 4:
3171  k = IA[i];
3172  j = JA[k];
3173  pA = val+k*25; // &val[k*jump];
3174  px0 = x+j*5; // &x[j*nb];
3175  py = py0;
3176  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3177 
3178  k ++;
3179  j = JA[k];
3180  pA = val+k*25; // &val[k*jump];
3181  px0 = x+j*5; // &x[j*nb];
3182  py = py0;
3183  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3184 
3185  k ++;
3186  j = JA[k];
3187  pA = val+k*25; // &val[k*jump];
3188  px0 = x+j*5; // &x[j*nb];
3189  py = py0;
3190  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3191 
3192  k ++;
3193  j = JA[k];
3194  pA = val+k*25; // &val[k*jump];
3195  px0 = x+j*5; // &x[j*nb];
3196  py = py0;
3197  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3198 
3199  break;
3200 
3201  case 5:
3202  k = IA[i];
3203  j = JA[k];
3204  pA = val+k*25; // &val[k*jump];
3205  px0 = x+j*5; // &x[j*nb];
3206  py = py0;
3207  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3208 
3209  k ++;
3210  j = JA[k];
3211  pA = val+k*25; // &val[k*jump];
3212  px0 = x+j*5; // &x[j*nb];
3213  py = py0;
3214  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3215 
3216  k ++;
3217  j = JA[k];
3218  pA = val+k*25; // &val[k*jump];
3219  px0 = x+j*5; // &x[j*nb];
3220  py = py0;
3221  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3222 
3223  k ++;
3224  j = JA[k];
3225  pA = val+k*25; // &val[k*jump];
3226  px0 = x+j*5; // &x[j*nb];
3227  py = py0;
3228  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3229 
3230  k ++;
3231  j = JA[k];
3232  pA = val+k*25; // &val[k*jump];
3233  px0 = x+j*5; // &x[j*nb];
3234  py = py0;
3235  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3236 
3237  break;
3238 
3239  case 6:
3240  k = IA[i];
3241  j = JA[k];
3242  pA = val+k*25; // &val[k*jump];
3243  px0 = x+j*5; // &x[j*nb];
3244  py = py0;
3245  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3246 
3247  k ++;
3248  j = JA[k];
3249  pA = val+k*25; // &val[k*jump];
3250  px0 = x+j*5; // &x[j*nb];
3251  py = py0;
3252  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3253 
3254  k ++;
3255  j = JA[k];
3256  pA = val+k*25; // &val[k*jump];
3257  px0 = x+j*5; // &x[j*nb];
3258  py = py0;
3259  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3260 
3261  k ++;
3262  j = JA[k];
3263  pA = val+k*25; // &val[k*jump];
3264  px0 = x+j*5; // &x[j*nb];
3265  py = py0;
3266  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3267 
3268  k ++;
3269  j = JA[k];
3270  pA = val+k*25; // &val[k*jump];
3271  px0 = x+j*5; // &x[j*nb];
3272  py = py0;
3273  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3274 
3275  k ++;
3276  j = JA[k];
3277  pA = val+k*25; // &val[k*jump];
3278  px0 = x+j*5; // &x[j*nb];
3279  py = py0;
3280  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3281 
3282  break;
3283 
3284  case 7:
3285  k = IA[i];
3286  j = JA[k];
3287  pA = val+k*25; // &val[k*jump];
3288  px0 = x+j*5; // &x[j*nb];
3289  py = py0;
3290  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3291 
3292  k ++;
3293  j = JA[k];
3294  pA = val+k*25; // &val[k*jump];
3295  px0 = x+j*5; // &x[j*nb];
3296  py = py0;
3297  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3298 
3299  k ++;
3300  j = JA[k];
3301  pA = val+k*25; // &val[k*jump];
3302  px0 = x+j*5; // &x[j*nb];
3303  py = py0;
3304  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3305 
3306  k ++;
3307  j = JA[k];
3308  pA = val+k*25; // &val[k*jump];
3309  px0 = x+j*5; // &x[j*nb];
3310  py = py0;
3311  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3312 
3313  k ++;
3314  j = JA[k];
3315  pA = val+k*25; // &val[k*jump];
3316  px0 = x+j*5; // &x[j*nb];
3317  py = py0;
3318  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3319 
3320  k ++;
3321  j = JA[k];
3322  pA = val+k*25; // &val[k*jump];
3323  px0 = x+j*5; // &x[j*nb];
3324  py = py0;
3325  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3326 
3327  k ++;
3328  j = JA[k];
3329  pA = val+k*25; // &val[k*jump];
3330  px0 = x+j*5; // &x[j*nb];
3331  py = py0;
3332  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3333 
3334  break;
3335 
3336  default:
3337  for (k = IA[i]; k < IA[i+1]; ++k)
3338  {
3339  j = JA[k];
3340  pA = val+k*25; // &val[k*jump];
3341  px0 = x+j*5; // &x[j*nb];
3342  py = py0;
3343  fasp_blas_smat_ypAx_nc5( pA, px0, py );
3344  }
3345  break;
3346  }
3347  }
3348  }
3349 #endif
3350  }
3351  else {
3352  for (i = 0; i < ROW; ++i) {
3353  py0 = &y[i*5];
3354  num_nnz_row = IA[i+1] - IA[i];
3355  switch(num_nnz_row) {
3356 
3357  case 3:
3358  k = IA[i];
3359  j = JA[k];
3360  px0 = x+j*5; // &x[j*nb];
3361  py = py0;
3362  py[0] += px0[0];
3363  py[1] += px0[1];
3364  py[2] += px0[2];
3365  py[3] += px0[3];
3366  py[4] += px0[4];
3367 
3368  k ++;
3369  j = JA[k];
3370  px0 = x+j*5; // &x[j*nb];
3371  py = py0;
3372  py[0] += px0[0];
3373  py[1] += px0[1];
3374  py[2] += px0[2];
3375  py[3] += px0[3];
3376  py[4] += px0[4];
3377 
3378  k ++;
3379  j = JA[k];
3380  px0 = x+j*5; // &x[j*nb];
3381  py = py0;
3382  py[0] += px0[0];
3383  py[1] += px0[1];
3384  py[2] += px0[2];
3385  py[3] += px0[3];
3386  py[4] += px0[4];
3387 
3388  break;
3389 
3390  case 4:
3391  k = IA[i];
3392  j = JA[k];
3393  px0 = x+j*5; // &x[j*nb];
3394  py = py0;
3395  py[0] += px0[0];
3396  py[1] += px0[1];
3397  py[2] += px0[2];
3398  py[3] += px0[3];
3399  py[4] += px0[4];
3400 
3401  k ++;
3402  j = JA[k];
3403  px0 = x+j*5; // &x[j*nb];
3404  py = py0;
3405  py[0] += px0[0];
3406  py[1] += px0[1];
3407  py[2] += px0[2];
3408  py[3] += px0[3];
3409  py[4] += px0[4];
3410 
3411  k ++;
3412  j = JA[k];
3413  px0 = x+j*5; // &x[j*nb];
3414  py = py0;
3415  py[0] += px0[0];
3416  py[1] += px0[1];
3417  py[2] += px0[2];
3418  py[3] += px0[3];
3419  py[4] += px0[4];
3420 
3421  k ++;
3422  j = JA[k];
3423  px0 = x+j*5; // &x[j*nb];
3424  py = py0;
3425  py[0] += px0[0];
3426  py[1] += px0[1];
3427  py[2] += px0[2];
3428  py[3] += px0[3];
3429  py[4] += px0[4];
3430 
3431  break;
3432 
3433  case 5:
3434  k = IA[i];
3435  j = JA[k];
3436  px0 = x+j*5; // &x[j*nb];
3437  py = py0;
3438  py[0] += px0[0];
3439  py[1] += px0[1];
3440  py[2] += px0[2];
3441  py[3] += px0[3];
3442  py[4] += px0[4];
3443 
3444  k ++;
3445  j = JA[k];
3446  px0 = x+j*5; // &x[j*nb];
3447  py = py0;
3448  py[0] += px0[0];
3449  py[1] += px0[1];
3450  py[2] += px0[2];
3451  py[3] += px0[3];
3452  py[4] += px0[4];
3453 
3454  k ++;
3455  j = JA[k];
3456  px0 = x+j*5; // &x[j*nb];
3457  py = py0;
3458  py[0] += px0[0];
3459  py[1] += px0[1];
3460  py[2] += px0[2];
3461  py[3] += px0[3];
3462  py[4] += px0[4];
3463 
3464  k ++;
3465  j = JA[k];
3466  px0 = x+j*5; // &x[j*nb];
3467  py = py0;
3468  py[0] += px0[0];
3469  py[1] += px0[1];
3470  py[2] += px0[2];
3471  py[3] += px0[3];
3472  py[4] += px0[4];
3473 
3474  k ++;
3475  j = JA[k];
3476  px0 = x+j*5; // &x[j*nb];
3477  py = py0;
3478  py[0] += px0[0];
3479  py[1] += px0[1];
3480  py[2] += px0[2];
3481  py[3] += px0[3];
3482  py[4] += px0[4];
3483 
3484  break;
3485 
3486  case 6:
3487  k = IA[i];
3488  j = JA[k];
3489  px0 = x+j*5; // &x[j*nb];
3490  py = py0;
3491  py[0] += px0[0];
3492  py[1] += px0[1];
3493  py[2] += px0[2];
3494  py[3] += px0[3];
3495  py[4] += px0[4];
3496 
3497  k ++;
3498  j = JA[k];
3499  px0 = x+j*5; // &x[j*nb];
3500  py = py0;
3501  py[0] += px0[0];
3502  py[1] += px0[1];
3503  py[2] += px0[2];
3504  py[3] += px0[3];
3505  py[4] += px0[4];
3506 
3507  k ++;
3508  j = JA[k];
3509  px0 = x+j*5; // &x[j*nb];
3510  py = py0;
3511  py[0] += px0[0];
3512  py[1] += px0[1];
3513  py[2] += px0[2];
3514  py[3] += px0[3];
3515  py[4] += px0[4];
3516 
3517  k ++;
3518  j = JA[k];
3519  px0 = x+j*5; // &x[j*nb];
3520  py = py0;
3521  py[0] += px0[0];
3522  py[1] += px0[1];
3523  py[2] += px0[2];
3524  py[3] += px0[3];
3525  py[4] += px0[4];
3526 
3527  k ++;
3528  j = JA[k];
3529  px0 = x+j*5; // &x[j*nb];
3530  py = py0;
3531  py[0] += px0[0];
3532  py[1] += px0[1];
3533  py[2] += px0[2];
3534  py[3] += px0[3];
3535  py[4] += px0[4];
3536 
3537  k ++;
3538  j = JA[k];
3539  px0 = x+j*5; // &x[j*nb];
3540  py = py0;
3541  py[0] += px0[0];
3542  py[1] += px0[1];
3543  py[2] += px0[2];
3544  py[3] += px0[3];
3545  py[4] += px0[4];
3546 
3547  break;
3548 
3549  case 7:
3550  k = IA[i];
3551  j = JA[k];
3552  px0 = x+j*5; // &x[j*nb];
3553  py = py0;
3554  py[0] += px0[0];
3555  py[1] += px0[1];
3556  py[2] += px0[2];
3557  py[3] += px0[3];
3558  py[4] += px0[4];
3559 
3560  k ++;
3561  j = JA[k];
3562  px0 = x+j*5; // &x[j*nb];
3563  py = py0;
3564  py[0] += px0[0];
3565  py[1] += px0[1];
3566  py[2] += px0[2];
3567  py[3] += px0[3];
3568  py[4] += px0[4];
3569 
3570  k ++;
3571  j = JA[k];
3572  px0 = x+j*5; // &x[j*nb];
3573  py = py0;
3574  py[0] += px0[0];
3575  py[1] += px0[1];
3576  py[2] += px0[2];
3577  py[3] += px0[3];
3578  py[4] += px0[4];
3579 
3580  k ++;
3581  j = JA[k];
3582  px0 = x+j*5; // &x[j*nb];
3583  py = py0;
3584  py[0] += px0[0];
3585  py[1] += px0[1];
3586  py[2] += px0[2];
3587  py[3] += px0[3];
3588  py[4] += px0[4];
3589 
3590  k ++;
3591  j = JA[k];
3592  px0 = x+j*5; // &x[j*nb];
3593  py = py0;
3594  py[0] += px0[0];
3595  py[1] += px0[1];
3596  py[2] += px0[2];
3597  py[3] += px0[3];
3598  py[4] += px0[4];
3599 
3600  k ++;
3601  j = JA[k];
3602  px0 = x+j*5; // &x[j*nb];
3603  py = py0;
3604  py[0] += px0[0];
3605  py[1] += px0[1];
3606  py[2] += px0[2];
3607  py[3] += px0[3];
3608  py[4] += px0[4];
3609 
3610  k ++;
3611  j = JA[k];
3612  px0 = x+j*5; // &x[j*nb];
3613  py = py0;
3614  py[0] += px0[0];
3615  py[1] += px0[1];
3616  py[2] += px0[2];
3617  py[3] += px0[3];
3618  py[4] += px0[4];
3619 
3620  break;
3621 
3622  default:
3623  for (k = IA[i]; k < IA[i+1]; ++k) {
3624  j = JA[k];
3625  px0 = x+j*5; // &x[j*nb];
3626  py = py0;
3627  py[0] += px0[0];
3628  py[1] += px0[1];
3629  py[2] += px0[2];
3630  py[3] += px0[3];
3631  py[4] += px0[4];
3632  }
3633  break;
3634  }
3635  }
3636  }
3637  }
3638  break;
3639 
3640  case 7:
3641  {
3642  if (use_openmp) {
3643 #ifdef _OPENMP
3644 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
3645  {
3646  myid = omp_get_thread_num();
3647  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
3648  for (i=mybegin; i < myend; ++i) {
3649  py0 = &y[i*7];
3650  num_nnz_row = IA[i+1] - IA[i];
3651  switch(num_nnz_row) {
3652 
3653  case 3:
3654  k = IA[i];
3655  j = JA[k];
3656  pA = val+k*49; // &val[k*jump];
3657  px0 = x+j*7; // &x[j*nb];
3658  py = py0;
3659  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3660 
3661  k ++;
3662  j = JA[k];
3663  pA = val+k*49; // &val[k*jump];
3664  px0 = x+j*7; // &x[j*nb];
3665  py = py0;
3666  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3667 
3668  k ++;
3669  j = JA[k];
3670  pA = val+k*49; // &val[k*jump];
3671  px0 = x+j*7; // &x[j*nb];
3672  py = py0;
3673  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3674 
3675  break;
3676 
3677  case 4:
3678  k = IA[i];
3679  j = JA[k];
3680  pA = val+k*49; // &val[k*jump];
3681  px0 = x+j*7; // &x[j*nb];
3682  py = py0;
3683  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3684 
3685  k ++;
3686  j = JA[k];
3687  pA = val+k*49; // &val[k*jump];
3688  px0 = x+j*7; // &x[j*nb];
3689  py = py0;
3690  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3691 
3692  k ++;
3693  j = JA[k];
3694  pA = val+k*49; // &val[k*jump];
3695  px0 = x+j*7; // &x[j*nb];
3696  py = py0;
3697  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3698 
3699  k ++;
3700  j = JA[k];
3701  pA = val+k*49; // &val[k*jump];
3702  px0 = x+j*7; // &x[j*nb];
3703  py = py0;
3704  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3705 
3706  break;
3707 
3708  case 5:
3709  k = IA[i];
3710  j = JA[k];
3711  pA = val+k*49; // &val[k*jump];
3712  px0 = x+j*7; // &x[j*nb];
3713  py = py0;
3714  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3715 
3716  k ++;
3717  j = JA[k];
3718  pA = val+k*49; // &val[k*jump];
3719  px0 = x+j*7; // &x[j*nb];
3720  py = py0;
3721  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3722 
3723  k ++;
3724  j = JA[k];
3725  pA = val+k*49; // &val[k*jump];
3726  px0 = x+j*7; // &x[j*nb];
3727  py = py0;
3728  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3729 
3730  k ++;
3731  j = JA[k];
3732  pA = val+k*49; // &val[k*jump];
3733  px0 = x+j*7; // &x[j*nb];
3734  py = py0;
3735  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3736 
3737  k ++;
3738  j = JA[k];
3739  pA = val+k*49; // &val[k*jump];
3740  px0 = x+j*7; // &x[j*nb];
3741  py = py0;
3742  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3743 
3744  break;
3745 
3746  case 6:
3747  k = IA[i];
3748  j = JA[k];
3749  pA = val+k*49; // &val[k*jump];
3750  px0 = x+j*7; // &x[j*nb];
3751  py = py0;
3752  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3753 
3754  k ++;
3755  j = JA[k];
3756  pA = val+k*49; // &val[k*jump];
3757  px0 = x+j*7; // &x[j*nb];
3758  py = py0;
3759  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3760 
3761  k ++;
3762  j = JA[k];
3763  pA = val+k*49; // &val[k*jump];
3764  px0 = x+j*7; // &x[j*nb];
3765  py = py0;
3766  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3767 
3768  k ++;
3769  j = JA[k];
3770  pA = val+k*49; // &val[k*jump];
3771  px0 = x+j*7; // &x[j*nb];
3772  py = py0;
3773  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3774 
3775  k ++;
3776  j = JA[k];
3777  pA = val+k*49; // &val[k*jump];
3778  px0 = x+j*7; // &x[j*nb];
3779  py = py0;
3780  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3781 
3782  k ++;
3783  j = JA[k];
3784  pA = val+k*49; // &val[k*jump];
3785  px0 = x+j*7; // &x[j*nb];
3786  py = py0;
3787  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3788 
3789  break;
3790 
3791  case 7:
3792  k = IA[i];
3793  j = JA[k];
3794  pA = val+k*49; // &val[k*jump];
3795  px0 = x+j*7; // &x[j*nb];
3796  py = py0;
3797  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3798 
3799  k ++;
3800  j = JA[k];
3801  pA = val+k*49; // &val[k*jump];
3802  px0 = x+j*7; // &x[j*nb];
3803  py = py0;
3804  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3805 
3806  k ++;
3807  j = JA[k];
3808  pA = val+k*49; // &val[k*jump];
3809  px0 = x+j*7; // &x[j*nb];
3810  py = py0;
3811  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3812 
3813  k ++;
3814  j = JA[k];
3815  pA = val+k*49; // &val[k*jump];
3816  px0 = x+j*7; // &x[j*nb];
3817  py = py0;
3818  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3819 
3820  k ++;
3821  j = JA[k];
3822  pA = val+k*49; // &val[k*jump];
3823  px0 = x+j*7; // &x[j*nb];
3824  py = py0;
3825  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3826 
3827  k ++;
3828  j = JA[k];
3829  pA = val+k*49; // &val[k*jump];
3830  px0 = x+j*7; // &x[j*nb];
3831  py = py0;
3832  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3833 
3834  k ++;
3835  j = JA[k];
3836  pA = val+k*49; // &val[k*jump];
3837  px0 = x+j*7; // &x[j*nb];
3838  py = py0;
3839  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3840 
3841  break;
3842 
3843  default:
3844  for (k = IA[i]; k < IA[i+1]; ++k) {
3845  j = JA[k];
3846  pA = val+k*49; // &val[k*jump];
3847  px0 = x+j*7; // &x[j*nb];
3848  py = py0;
3849  fasp_blas_smat_ypAx_nc7( pA, px0, py );
3850  }
3851  break;
3852  }
3853  }
3854  }
3855 #endif
3856  }
3857  else {
3858  for (i = 0; i < ROW; ++i) {
3859  py0 = &y[i*7];
3860  num_nnz_row = IA[i+1] - IA[i];
3861  switch(num_nnz_row) {
3862 
3863  case 3:
3864  k = IA[i];
3865  j = JA[k];
3866  px0 = x+j*7; // &x[j*nb];
3867  py = py0;
3868  py[0] += px0[0];
3869  py[1] += px0[1];
3870  py[2] += px0[2];
3871  py[3] += px0[3];
3872  py[4] += px0[4];
3873  py[5] += px0[5];
3874  py[6] += px0[6];
3875 
3876  k ++;
3877  j = JA[k];
3878  px0 = x+j*7; // &x[j*nb];
3879  py = py0;
3880  py[0] += px0[0];
3881  py[1] += px0[1];
3882  py[2] += px0[2];
3883  py[3] += px0[3];
3884  py[4] += px0[4];
3885  py[5] += px0[5];
3886  py[6] += px0[6];
3887 
3888  k ++;
3889  j = JA[k];
3890  px0 = x+j*7; // &x[j*nb];
3891  py = py0;
3892  py[0] += px0[0];
3893  py[1] += px0[1];
3894  py[2] += px0[2];
3895  py[3] += px0[3];
3896  py[4] += px0[4];
3897  py[5] += px0[5];
3898  py[6] += px0[6];
3899 
3900  break;
3901 
3902  case 4:
3903  k = IA[i];
3904  j = JA[k];
3905  px0 = x+j*7; // &x[j*nb];
3906  py = py0;
3907  py[0] += px0[0];
3908  py[1] += px0[1];
3909  py[2] += px0[2];
3910  py[3] += px0[3];
3911  py[4] += px0[4];
3912  py[5] += px0[5];
3913  py[6] += px0[6];
3914 
3915  k ++;
3916  j = JA[k];
3917  px0 = x+j*7; // &x[j*nb];
3918  py = py0;
3919  py[0] += px0[0];
3920  py[1] += px0[1];
3921  py[2] += px0[2];
3922  py[3] += px0[3];
3923  py[4] += px0[4];
3924  py[5] += px0[5];
3925  py[6] += px0[6];
3926 
3927  k ++;
3928  j = JA[k];
3929  px0 = x+j*7; // &x[j*nb];
3930  py = py0;
3931  py[0] += px0[0];
3932  py[1] += px0[1];
3933  py[2] += px0[2];
3934  py[3] += px0[3];
3935  py[4] += px0[4];
3936  py[5] += px0[5];
3937  py[6] += px0[6];
3938 
3939  k ++;
3940  j = JA[k];
3941  px0 = x+j*7; // &x[j*nb];
3942  py = py0;
3943  py[0] += px0[0];
3944  py[1] += px0[1];
3945  py[2] += px0[2];
3946  py[3] += px0[3];
3947  py[4] += px0[4];
3948  py[5] += px0[5];
3949  py[6] += px0[6];
3950 
3951  break;
3952 
3953  case 5:
3954  k = IA[i];
3955  j = JA[k];
3956  px0 = x+j*7; // &x[j*nb];
3957  py = py0;
3958  py[0] += px0[0];
3959  py[1] += px0[1];
3960  py[2] += px0[2];
3961  py[3] += px0[3];
3962  py[4] += px0[4];
3963  py[5] += px0[5];
3964  py[6] += px0[6];
3965 
3966  k ++;
3967  j = JA[k];
3968  px0 = x+j*7; // &x[j*nb];
3969  py = py0;
3970  py[0] += px0[0];
3971  py[1] += px0[1];
3972  py[2] += px0[2];
3973  py[3] += px0[3];
3974  py[4] += px0[4];
3975  py[5] += px0[5];
3976  py[6] += px0[6];
3977 
3978  k ++;
3979  j = JA[k];
3980  px0 = x+j*7; // &x[j*nb];
3981  py = py0;
3982  py[0] += px0[0];
3983  py[1] += px0[1];
3984  py[2] += px0[2];
3985  py[3] += px0[3];
3986  py[4] += px0[4];
3987  py[5] += px0[5];
3988  py[6] += px0[6];
3989 
3990  k ++;
3991  j = JA[k];
3992  px0 = x+j*7; // &x[j*nb];
3993  py = py0;
3994  py[0] += px0[0];
3995  py[1] += px0[1];
3996  py[2] += px0[2];
3997  py[3] += px0[3];
3998  py[4] += px0[4];
3999  py[5] += px0[5];
4000  py[6] += px0[6];
4001 
4002  k ++;
4003  j = JA[k];
4004  px0 = x+j*7; // &x[j*nb];
4005  py = py0;
4006  py[0] += px0[0];
4007  py[1] += px0[1];
4008  py[2] += px0[2];
4009  py[3] += px0[3];
4010  py[4] += px0[4];
4011  py[5] += px0[5];
4012  py[6] += px0[6];
4013 
4014  break;
4015 
4016  case 6:
4017  k = IA[i];
4018  j = JA[k];
4019  px0 = x+j*7; // &x[j*nb];
4020  py = py0;
4021  py[0] += px0[0];
4022  py[1] += px0[1];
4023  py[2] += px0[2];
4024  py[3] += px0[3];
4025  py[4] += px0[4];
4026  py[5] += px0[5];
4027  py[6] += px0[6];
4028 
4029  k ++;
4030  j = JA[k];
4031  px0 = x+j*7; // &x[j*nb];
4032  py = py0;
4033  py[0] += px0[0];
4034  py[1] += px0[1];
4035  py[2] += px0[2];
4036  py[3] += px0[3];
4037  py[4] += px0[4];
4038  py[5] += px0[5];
4039  py[6] += px0[6];
4040 
4041  k ++;
4042  j = JA[k];
4043  px0 = x+j*7; // &x[j*nb];
4044  py = py0;
4045  py[0] += px0[0];
4046  py[1] += px0[1];
4047  py[2] += px0[2];
4048  py[3] += px0[3];
4049  py[4] += px0[4];
4050  py[5] += px0[5];
4051  py[6] += px0[6];
4052 
4053  k ++;
4054  j = JA[k];
4055  px0 = x+j*7; // &x[j*nb];
4056  py = py0;
4057  py[0] += px0[0];
4058  py[1] += px0[1];
4059  py[2] += px0[2];
4060  py[3] += px0[3];
4061  py[4] += px0[4];
4062  py[5] += px0[5];
4063  py[6] += px0[6];
4064 
4065  k ++;
4066  j = JA[k];
4067  px0 = x+j*7; // &x[j*nb];
4068  py = py0;
4069  py[0] += px0[0];
4070  py[1] += px0[1];
4071  py[2] += px0[2];
4072  py[3] += px0[3];
4073  py[4] += px0[4];
4074  py[5] += px0[5];
4075  py[6] += px0[6];
4076 
4077  k ++;
4078  j = JA[k];
4079  px0 = x+j*7; // &x[j*nb];
4080  py = py0;
4081  py[0] += px0[0];
4082  py[1] += px0[1];
4083  py[2] += px0[2];
4084  py[3] += px0[3];
4085  py[4] += px0[4];
4086  py[5] += px0[5];
4087  py[6] += px0[6];
4088 
4089  break;
4090 
4091  case 7:
4092  k = IA[i];
4093  j = JA[k];
4094  px0 = x+j*7; // &x[j*nb];
4095  py = py0;
4096  py[0] += px0[0];
4097  py[1] += px0[1];
4098  py[2] += px0[2];
4099  py[3] += px0[3];
4100  py[4] += px0[4];
4101  py[5] += px0[5];
4102  py[6] += px0[6];
4103 
4104  k ++;
4105  j = JA[k];
4106  px0 = x+j*7; // &x[j*nb];
4107  py = py0;
4108  py[0] += px0[0];
4109  py[1] += px0[1];
4110  py[2] += px0[2];
4111  py[3] += px0[3];
4112  py[4] += px0[4];
4113  py[5] += px0[5];
4114  py[6] += px0[6];
4115 
4116  k ++;
4117  j = JA[k];
4118  px0 = x+j*7; // &x[j*nb];
4119  py = py0;
4120  py[0] += px0[0];
4121  py[1] += px0[1];
4122  py[2] += px0[2];
4123  py[3] += px0[3];
4124  py[4] += px0[4];
4125  py[5] += px0[5];
4126  py[6] += px0[6];
4127 
4128  k ++;
4129  j = JA[k];
4130  px0 = x+j*7; // &x[j*nb];
4131  py = py0;
4132  py[0] += px0[0];
4133  py[1] += px0[1];
4134  py[2] += px0[2];
4135  py[3] += px0[3];
4136  py[4] += px0[4];
4137  py[5] += px0[5];
4138  py[6] += px0[6];
4139 
4140  k ++;
4141  j = JA[k];
4142  px0 = x+j*7; // &x[j*nb];
4143  py = py0;
4144  py[0] += px0[0];
4145  py[1] += px0[1];
4146  py[2] += px0[2];
4147  py[3] += px0[3];
4148  py[4] += px0[4];
4149  py[5] += px0[5];
4150  py[6] += px0[6];
4151 
4152  k ++;
4153  j = JA[k];
4154  px0 = x+j*7; // &x[j*nb];
4155  py = py0;
4156  py[0] += px0[0];
4157  py[1] += px0[1];
4158  py[2] += px0[2];
4159  py[3] += px0[3];
4160  py[4] += px0[4];
4161  py[5] += px0[5];
4162  py[6] += px0[6];
4163 
4164  k ++;
4165  j = JA[k];
4166  px0 = x+j*7; // &x[j*nb];
4167  py = py0;
4168  py[0] += px0[0];
4169  py[1] += px0[1];
4170  py[2] += px0[2];
4171  py[3] += px0[3];
4172  py[4] += px0[4];
4173  py[5] += px0[5];
4174  py[6] += px0[6];
4175 
4176  break;
4177 
4178  default:
4179  for (k = IA[i]; k < IA[i+1]; ++k) {
4180  j = JA[k];
4181  px0 = x+j*7; // &x[j*nb];
4182  py = py0;
4183  py[0] += px0[0];
4184  py[1] += px0[1];
4185  py[2] += px0[2];
4186  py[3] += px0[3];
4187  py[4] += px0[4];
4188  py[5] += px0[5];
4189  py[6] += px0[6];
4190  }
4191  break;
4192  }
4193  }
4194  }
4195  }
4196  break;
4197 
4198  default:
4199  {
4200  if (use_openmp) {
4201 #ifdef _OPENMP
4202 #pragma omp parallel private(myid, mybegin, myend, i, py0, num_nnz_row, k, j, pA, px0, py)
4203  {
4204  myid = omp_get_thread_num();
4205  FASP_GET_START_END(myid, nthreads, ROW, &mybegin, &myend);
4206  for (i=mybegin; i < myend; ++i) {
4207  py0 = &y[i*nb];
4208  num_nnz_row = IA[i+1] - IA[i];
4209  switch(num_nnz_row) {
4210 
4211  case 3:
4212  k = IA[i];
4213  j = JA[k];
4214  px0 = x+j*nb; // &x[j*nb];
4215  py = py0;
4216  fasp_blas_array_axpy (nb, 1.0, px0, py);
4217 
4218  k ++;
4219  j = JA[k];
4220  px0 = x+j*nb; // &x[j*nb];
4221  py = py0;
4222  fasp_blas_array_axpy (nb, 1.0, px0, py);
4223 
4224  k ++;
4225  j = JA[k];
4226  px0 = x+j*nb; // &x[j*nb];
4227  py = py0;
4228  fasp_blas_array_axpy (nb, 1.0, px0, py);
4229 
4230  break;
4231 
4232  case 4:
4233  k = IA[i];
4234  j = JA[k];
4235  px0 = x+j*nb; // &x[j*nb];
4236  py = py0;
4237  fasp_blas_array_axpy (nb, 1.0, px0, py);
4238 
4239  k ++;
4240  j = JA[k];
4241  px0 = x+j*nb; // &x[j*nb];
4242  py = py0;
4243  fasp_blas_array_axpy (nb, 1.0, px0, py);
4244 
4245  k ++;
4246  j = JA[k];
4247  px0 = x+j*nb; // &x[j*nb];
4248  py = py0;
4249  fasp_blas_array_axpy (nb, 1.0, px0, py);
4250 
4251  k ++;
4252  j = JA[k];
4253  px0 = x+j*nb; // &x[j*nb];
4254  py = py0;
4255  fasp_blas_array_axpy (nb, 1.0, px0, py);
4256 
4257  break;
4258 
4259  case 5:
4260  k = IA[i];
4261  j = JA[k];
4262  px0 = x+j*nb; // &x[j*nb];
4263  py = py0;
4264  fasp_blas_array_axpy (nb, 1.0, px0, py);
4265 
4266  k ++;
4267  j = JA[k];
4268  px0 = x+j*nb; // &x[j*nb];
4269  py = py0;
4270  fasp_blas_array_axpy (nb, 1.0, px0, py);
4271 
4272  k ++;
4273  j = JA[k];
4274  px0 = x+j*nb; // &x[j*nb];
4275  py = py0;
4276  fasp_blas_array_axpy (nb, 1.0, px0, py);
4277 
4278  k ++;
4279  j = JA[k];
4280  px0 = x+j*nb; // &x[j*nb];
4281  py = py0;
4282  fasp_blas_array_axpy (nb, 1.0, px0, py);
4283 
4284  k ++;
4285  j = JA[k];
4286  px0 = x+j*nb; // &x[j*nb];
4287  py = py0;
4288  fasp_blas_array_axpy (nb, 1.0, px0, py);
4289 
4290  break;
4291 
4292  case 6:
4293  k = IA[i];
4294  j = JA[k];
4295  px0 = x+j*nb; // &x[j*nb];
4296  py = py0;
4297  fasp_blas_array_axpy (nb, 1.0, px0, py);
4298 
4299  k ++;
4300  j = JA[k];
4301  px0 = x+j*nb; // &x[j*nb];
4302  py = py0;
4303  fasp_blas_array_axpy (nb, 1.0, px0, py);
4304 
4305  k ++;
4306  j = JA[k];
4307  px0 = x+j*nb; // &x[j*nb];
4308  py = py0;
4309  fasp_blas_array_axpy (nb, 1.0, px0, py);
4310 
4311  k ++;
4312  j = JA[k];
4313  px0 = x+j*nb; // &x[j*nb];
4314  py = py0;
4315  fasp_blas_array_axpy (nb, 1.0, px0, py);
4316 
4317  k ++;
4318  j = JA[k];
4319  px0 = x+j*nb; // &x[j*nb];
4320  py = py0;
4321  fasp_blas_array_axpy (nb, 1.0, px0, py);
4322 
4323  k ++;
4324  j = JA[k];
4325  px0 = x+j*nb; // &x[j*nb];
4326  py = py0;
4327  fasp_blas_array_axpy (nb, 1.0, px0, py);
4328 
4329  break;
4330 
4331  case 7:
4332  k = IA[i];
4333  j = JA[k];
4334  px0 = x+j*nb; // &x[j*nb];
4335  py = py0;
4336  fasp_blas_array_axpy (nb, 1.0, px0, py);
4337 
4338  k ++;
4339  j = JA[k];
4340  px0 = x+j*nb; // &x[j*nb];
4341  py = py0;
4342  fasp_blas_array_axpy (nb, 1.0, px0, py);
4343 
4344  k ++;
4345  j = JA[k];
4346  px0 = x+j*nb; // &x[j*nb];
4347  py = py0;
4348  fasp_blas_array_axpy (nb, 1.0, px0, py);
4349 
4350  k ++;
4351  j = JA[k];
4352  px0 = x+j*nb; // &x[j*nb];
4353  py = py0;
4354  fasp_blas_array_axpy (nb, 1.0, px0, py);
4355 
4356  k ++;
4357  j = JA[k];
4358  px0 = x+j*nb; // &x[j*nb];
4359  py = py0;
4360  fasp_blas_array_axpy (nb, 1.0, px0, py);
4361 
4362  k ++;
4363  j = JA[k];
4364  px0 = x+j*nb; // &x[j*nb];
4365  py = py0;
4366  fasp_blas_array_axpy (nb, 1.0, px0, py);
4367 
4368  k ++;
4369  j = JA[k];
4370  px0 = x+j*nb; // &x[j*nb];
4371  py = py0;
4372  fasp_blas_array_axpy (nb, 1.0, px0, py);
4373 
4374  break;
4375 
4376  default:
4377  for (k = IA[i]; k < IA[i+1]; ++k) {
4378  j = JA[k];
4379  px0 = x+j*nb; // &x[j*nb];
4380  py = py0;
4381  fasp_blas_array_axpy (nb, 1.0, px0, py);
4382  }
4383  break;
4384  }
4385  }
4386  }
4387 #endif
4388  }
4389  else {
4390  for (i = 0; i < ROW; ++i) {
4391  py0 = &y[i*nb];
4392  num_nnz_row = IA[i+1] - IA[i];
4393  switch(num_nnz_row) {
4394 
4395  case 3:
4396  k = IA[i];
4397  j = JA[k];
4398  px0 = x+j*nb; // &x[j*nb];
4399  py = py0;
4400  fasp_blas_array_axpy (nb, 1.0, px0, py);
4401 
4402  k ++;
4403  j = JA[k];
4404  px0 = x+j*nb; // &x[j*nb];
4405  py = py0;
4406  fasp_blas_array_axpy (nb, 1.0, px0, py);
4407 
4408  k ++;
4409  j = JA[k];
4410  px0 = x+j*nb; // &x[j*nb];
4411  py = py0;
4412  fasp_blas_array_axpy (nb, 1.0, px0, py);
4413 
4414  break;
4415 
4416  case 4:
4417  k = IA[i];
4418  j = JA[k];
4419  px0 = x+j*nb; // &x[j*nb];
4420  py = py0;
4421  fasp_blas_array_axpy (nb, 1.0, px0, py);
4422 
4423  k ++;
4424  j = JA[k];
4425  px0 = x+j*nb; // &x[j*nb];
4426  py = py0;
4427  fasp_blas_array_axpy (nb, 1.0, px0, py);
4428 
4429  k ++;
4430  j = JA[k];
4431  px0 = x+j*nb; // &x[j*nb];
4432  py = py0;
4433  fasp_blas_array_axpy (nb, 1.0, px0, py);
4434 
4435  k ++;
4436  j = JA[k];
4437  px0 = x+j*nb; // &x[j*nb];
4438  py = py0;
4439  fasp_blas_array_axpy (nb, 1.0, px0, py);
4440 
4441  break;
4442 
4443  case 5:
4444  k = IA[i];
4445  j = JA[k];
4446  px0 = x+j*nb; // &x[j*nb];
4447  py = py0;
4448  fasp_blas_array_axpy (nb, 1.0, px0, py);
4449 
4450  k ++;
4451  j = JA[k];
4452  px0 = x+j*nb; // &x[j*nb];
4453  py = py0;
4454  fasp_blas_array_axpy (nb, 1.0, px0, py);
4455 
4456  k ++;
4457  j = JA[k];
4458  px0 = x+j*nb; // &x[j*nb];
4459  py = py0;
4460  fasp_blas_array_axpy (nb, 1.0, px0, py);
4461 
4462  k ++;
4463  j = JA[k];
4464  px0 = x+j*nb; // &x[j*nb];
4465  py = py0;
4466  fasp_blas_array_axpy (nb, 1.0, px0, py);
4467 
4468  k ++;
4469  j = JA[k];
4470  px0 = x+j*nb; // &x[j*nb];
4471  py = py0;
4472  fasp_blas_array_axpy (nb, 1.0, px0, py);
4473 
4474  break;
4475 
4476  case 6:
4477  k = IA[i];
4478  j = JA[k];
4479  px0 = x+j*nb; // &x[j*nb];
4480  py = py0;
4481  fasp_blas_array_axpy (nb, 1.0, px0, py);
4482 
4483  k ++;
4484  j = JA[k];
4485  px0 = x+j*nb; // &x[j*nb];
4486  py = py0;
4487  fasp_blas_array_axpy (nb, 1.0, px0, py);
4488 
4489  k ++;
4490  j = JA[k];
4491  px0 = x+j*nb; // &x[j*nb];
4492  py = py0;
4493  fasp_blas_array_axpy (nb, 1.0, px0, py);
4494 
4495  k ++;
4496  j = JA[k];
4497  px0 = x+j*nb; // &x[j*nb];
4498  py = py0;
4499  fasp_blas_array_axpy (nb, 1.0, px0, py);
4500 
4501  k ++;
4502  j = JA[k];
4503  px0 = x+j*nb; // &x[j*nb];
4504  py = py0;
4505  fasp_blas_array_axpy (nb, 1.0, px0, py);
4506 
4507  k ++;
4508  j = JA[k];
4509  px0 = x+j*nb; // &x[j*nb];
4510  py = py0;
4511  fasp_blas_array_axpy (nb, 1.0, px0, py);
4512 
4513  break;
4514 
4515  case 7:
4516  k = IA[i];
4517  j = JA[k];
4518  px0 = x+j*nb; // &x[j*nb];
4519  py = py0;
4520  fasp_blas_array_axpy (nb, 1.0, px0, py);
4521 
4522  k ++;
4523  j = JA[k];
4524  px0 = x+j*nb; // &x[j*nb];
4525  py = py0;
4526  fasp_blas_array_axpy (nb, 1.0, px0, py);
4527 
4528  k ++;
4529  j = JA[k];
4530  px0 = x+j*nb; // &x[j*nb];
4531  py = py0;
4532  fasp_blas_array_axpy (nb, 1.0, px0, py);
4533 
4534  k ++;
4535  j = JA[k];
4536  px0 = x+j*nb; // &x[j*nb];
4537  py = py0;
4538  fasp_blas_array_axpy (nb, 1.0, px0, py);
4539 
4540  k ++;
4541  j = JA[k];
4542  px0 = x+j*nb; // &x[j*nb];
4543  py = py0;
4544  fasp_blas_array_axpy (nb, 1.0, px0, py);
4545 
4546  k ++;
4547  j = JA[k];
4548  px0 = x+j*nb; // &x[j*nb];
4549  py = py0;
4550  fasp_blas_array_axpy (nb, 1.0, px0, py);
4551 
4552  k ++;
4553  j = JA[k];
4554  px0 = x+j*nb; // &x[j*nb];
4555  py = py0;
4556  fasp_blas_array_axpy (nb, 1.0, px0, py);
4557 
4558  break;
4559 
4560  default:
4561  for (k = IA[i]; k < IA[i+1]; ++k) {
4562  j = JA[k];
4563  px0 = x+j*nb; // &x[j*nb];
4564  py = py0;
4565  fasp_blas_array_axpy (nb, 1.0, px0, py);
4566  }
4567  break;
4568  }
4569  }
4570  }
4571  }
4572  break;
4573  }
4574 }
4575 
4576 
4592  dBSRmat *B,
4593  dBSRmat *C)
4594 {
4595 
4596  INT i,j,k,l,count;
4597  INT *JD = (INT *)fasp_mem_calloc(B->COL,sizeof(INT));
4598 
4599  const INT nb = A->nb;
4600  const INT nb2 = nb*nb;
4601 
4602  // /TODO: check A and B see if there are compatible for multiplication!! -- Xiaozhe
4603 
4604  C->ROW = A->ROW;
4605  C->COL = B->COL;
4606  C->nb = A->nb;
4608 
4609  C->val = NULL;
4610  C->JA = NULL;
4611  C->IA = (INT*)fasp_mem_calloc(C->ROW+1,sizeof(INT));
4612 
4613  REAL *temp = (REAL *)fasp_mem_calloc(nb2, sizeof(REAL));
4614 
4615  for (i=0;i<B->COL;++i) JD[i]=-1;
4616 
4617  // step 1: Find first the structure IA of C
4618  for(i=0;i<C->ROW;++i) {
4619  count=0;
4620 
4621  for (k=A->IA[i];k<A->IA[i+1];++k) {
4622  for (j=B->IA[A->JA[k]];j<B->IA[A->JA[k]+1];++j) {
4623  for (l=0;l<count;l++) {
4624  if (JD[l]==B->JA[j]) break;
4625  }
4626 
4627  if (l==count) {
4628  JD[count]=B->JA[j];
4629  count++;
4630  }
4631  }
4632  }
4633  C->IA[i+1]=count;
4634  for (j=0;j<count;++j) {
4635  JD[j]=-1;
4636  }
4637  }
4638 
4639  for (i=0;i<C->ROW;++i) C->IA[i+1]+=C->IA[i];
4640 
4641  // step 2: Find the structure JA of C
4642  INT countJD;
4643 
4644  C->JA=(INT*)fasp_mem_calloc(C->IA[C->ROW],sizeof(INT));
4645 
4646  for (i=0;i<C->ROW;++i) {
4647  countJD=0;
4648  count=C->IA[i];
4649  for (k=A->IA[i];k<A->IA[i+1];++k) {
4650  for (j=B->IA[A->JA[k]];j<B->IA[A->JA[k]+1];++j) {
4651  for (l=0;l<countJD;l++) {
4652  if (JD[l]==B->JA[j]) break;
4653  }
4654 
4655  if (l==countJD) {
4656  C->JA[count]=B->JA[j];
4657  JD[countJD]=B->JA[j];
4658  count++;
4659  countJD++;
4660  }
4661  }
4662  }
4663 
4664  //for (j=0;j<countJD;++j) JD[j]=-1;
4665  fasp_iarray_set(countJD, JD, -1);
4666  }
4667 
4668  fasp_mem_free(JD);
4669 
4670  // step 3: Find the structure A of C
4671  C->val=(REAL*)fasp_mem_calloc((C->IA[C->ROW])*nb2,sizeof(REAL));
4672 
4673  for (i=0;i<C->ROW;++i) {
4674  for (j=C->IA[i];j<C->IA[i+1];++j) {
4675 
4676  fasp_array_set(nb2, C->val+(j*nb2), 0x0);
4677 
4678  for (k=A->IA[i];k<A->IA[i+1];++k) {
4679  for (l=B->IA[A->JA[k]];l<B->IA[A->JA[k]+1];l++) {
4680  if (B->JA[l]==C->JA[j]) {
4681  fasp_blas_smat_mul (A->val+(k*nb2), B->val+(l*nb2), temp, nb);
4682  fasp_blas_array_axpy (nb2, 1.0, temp, C->val+(j*nb2));
4683  } // end if
4684  } // end for l
4685  } // end for k
4686  } // end for j
4687  } // end for i
4688 
4689  C->NNZ = C->IA[C->ROW]-C->IA[0];
4690 
4691  fasp_mem_free(temp);
4692 
4693 }
4694 
4712  dBSRmat *A,
4713  dBSRmat *P,
4714  dBSRmat *B)
4715 {
4716  const INT row=R->ROW, col=P->COL,nb=A->nb, nb2=A->nb*A->nb;
4717  INT nB=A->NNZ;
4718  INT i,i1,j,jj,k,length;
4719  INT begin_row,end_row,begin_rowA,end_rowA,begin_rowR,end_rowR;
4720  INT istart,iistart,count;
4721 
4722  REAL *rj=R->val, *aj=A->val, *pj=P->val, *acj;
4723  INT *ir=R->IA, *ia=A->IA, *ip=P->IA, *iac;
4724  INT *jr=R->JA, *ja=A->JA, *jp=P->JA, *jac;
4725 
4726  INT *index=(INT *)fasp_mem_calloc(A->COL,sizeof(INT));
4727 
4728  REAL *smat_tmp=(REAL *)fasp_mem_calloc(nb2,sizeof(REAL));
4729 
4730  INT *iindex=(INT *)fasp_mem_calloc(col,sizeof(INT));
4731 
4732  for (i=0; i<A->COL; ++i) index[i] = -2;
4733 
4734  memcpy(iindex,index,col*sizeof(INT));
4735 
4736  jac=(INT*)fasp_mem_calloc(nB,sizeof(INT));
4737 
4738  iac=(INT*)fasp_mem_calloc(row+1,sizeof(INT));
4739 
4740  REAL *temp=(REAL*)fasp_mem_calloc(A->COL*nb2,sizeof(REAL));
4741 
4742  iac[0] = 0;
4743 
4744  // First loop: form sparsity pattern of R*A*P
4745  for (i=0; i < row; ++i) {
4746  // reset istart and length at the beginning of each loop
4747  istart = -1; length = 0; i1 = i+1;
4748 
4749  // go across the rows in R
4750  begin_rowR=ir[i]; end_rowR=ir[i1];
4751  for (jj=begin_rowR; jj<end_rowR; ++jj) {
4752  j = jr[jj];
4753  // for each column in A
4754  begin_rowA=ia[j]; end_rowA=ia[j+1];
4755  for (k=begin_rowA; k<end_rowA; ++k) {
4756  if (index[ja[k]] == -2) {
4757  index[ja[k]] = istart;
4758  istart = ja[k];
4759  ++length;
4760  }
4761  }
4762  }
4763 
4764  // book-keeping [resetting length and setting iistart]
4765  count = length; iistart = -1; length = 0;
4766 
4767  // use each column that would have resulted from R*A
4768  for (j=0; j < count; ++j) {
4769  jj = istart;
4770  istart = index[istart];
4771  index[jj] = -2;
4772 
4773  // go across the row of P
4774  begin_row=ip[jj]; end_row=ip[jj+1];
4775  for (k=begin_row; k<end_row; ++k) {
4776  // pull out the appropriate columns of P
4777  if (iindex[jp[k]] == -2){
4778  iindex[jp[k]] = iistart;
4779  iistart = jp[k];
4780  ++length;
4781  }
4782  } // end for k
4783  } // end for j
4784 
4785  // set B->IA
4786  iac[i1]=iac[i]+length;
4787 
4788  if (iac[i1]>nB) {
4789  nB=nB*2;
4790  jac=(INT*)fasp_mem_realloc(jac, nB*sizeof(INT));
4791  }
4792 
4793  // put the correct columns of p into the column list of the products
4794  begin_row=iac[i]; end_row=iac[i1];
4795  for (j=begin_row; j<end_row; ++j) {
4796  // put the value in B->JA
4797  jac[j] = iistart;
4798  // set istart to the next value
4799  iistart = iindex[iistart];
4800  // set the iindex spot to 0
4801  iindex[jac[j]] = -2;
4802  } // end j
4803  } // end i: First loop
4804 
4805  jac=(INT*)fasp_mem_realloc(jac,(iac[row])*sizeof(INT));
4806 
4807  acj=(REAL*)fasp_mem_calloc(iac[row]*nb2,sizeof(REAL));
4808 
4809  INT *BTindex=(INT*)fasp_mem_calloc(col,sizeof(INT));
4810 
4811  // Second loop: compute entries of R*A*P
4812  for (i=0; i<row; ++i) {
4813  i1 = i+1;
4814  // each col of B
4815  begin_row=iac[i]; end_row=iac[i1];
4816  for (j=begin_row; j<end_row; ++j) {
4817  BTindex[jac[j]]=j;
4818  }
4819  // reset istart and length at the beginning of each loop
4820  istart = -1; length = 0;
4821 
4822  // go across the rows in R
4823  begin_rowR=ir[i]; end_rowR=ir[i1];
4824  for ( jj=begin_rowR; jj<end_rowR; ++jj ) {
4825  j = jr[jj];
4826  // for each column in A
4827  begin_rowA=ia[j]; end_rowA=ia[j+1];
4828  for (k=begin_rowA; k<end_rowA; ++k) {
4829  if (index[ja[k]] == -2) {
4830  index[ja[k]] = istart;
4831  istart = ja[k];
4832  ++length;
4833  }
4834  fasp_blas_smat_mul(&rj[jj*nb2],&aj[k*nb2],smat_tmp,nb);
4835  //fasp_array_xpy(nb2,&temp[ja[k]*nb2], smat_tmp );
4836  fasp_blas_array_axpy (nb2, 1.0, smat_tmp, &temp[ja[k]*nb2]);
4837 
4838  //temp[ja[k]]+=rj[jj]*aj[k];
4839  // change to X = X+Y*Z
4840  }
4841  }
4842 
4843  // book-keeping [resetting length and setting iistart]
4844  // use each column that would have resulted from R*A
4845  for (j=0; j<length; ++j) {
4846  jj = istart;
4847  istart = index[istart];
4848  index[jj] = -2;
4849 
4850  // go across the row of P
4851  begin_row=ip[jj]; end_row=ip[jj+1];
4852  for (k=begin_row; k<end_row; ++k) {
4853  // pull out the appropriate columns of P
4854  //acj[BTindex[jp[k]]]+=temp[jj]*pj[k];
4855  fasp_blas_smat_mul(&temp[jj*nb2],&pj[k*nb2],smat_tmp,nb);
4856  //fasp_array_xpy(nb2,&acj[BTindex[jp[k]]*nb2], smat_tmp );
4857  fasp_blas_array_axpy (nb2, 1.0, smat_tmp, &acj[BTindex[jp[k]]*nb2]);
4858 
4859  // change to X = X+Y*Z
4860  }
4861  //temp[jj]=0.0; // change to X[nb,nb] = 0;
4862  fasp_array_set(nb2,&temp[jj*nb2],0.0);
4863  }
4864  } // end for i: Second loop
4865  // setup coarse matrix B
4866  B->ROW=row; B->COL=col;
4867  B->IA=iac; B->JA=jac; B->val=acj;
4868  B->NNZ=B->IA[B->ROW]-B->IA[0];
4869 
4870  B->nb=A->nb;
4872  fasp_mem_free(temp);
4873  fasp_mem_free(index);
4874  fasp_mem_free(iindex);
4875  fasp_mem_free(BTindex);
4876  fasp_mem_free(smat_tmp);
4877 }
4878 
4896  dBSRmat *A,
4897  dBSRmat *P,
4898  dBSRmat *B)
4899 {
4900  const INT row=R->ROW, col=P->COL,nb=A->nb, nb2=A->nb*A->nb;
4901  //INT nB=A->NNZ;
4902 
4903  REAL *rj=R->val, *aj=A->val, *pj=P->val, *acj;
4904  INT *ir=R->IA, *ia=A->IA, *ip=P->IA, *iac;
4905  INT *jr=R->JA, *ja=A->JA, *jp=P->JA, *jac;
4906 
4907  INT *Ps_marker = NULL;
4908  INT *As_marker = NULL;
4909 
4910 #ifdef _OPENMP
4911  INT *P_marker = NULL;
4912  INT *A_marker = NULL;
4913  REAL *smat_tmp = NULL;
4914 #endif
4915 
4916  INT i, i1, i2, i3, jj1, jj2, jj3;
4917  INT counter, jj_row_begining;
4918 
4919  INT nthreads = 1;
4920 
4921 #ifdef _OPENMP
4922  INT myid, mybegin, myend, Ctemp;
4923  nthreads = FASP_GET_NUM_THREADS();
4924 #endif
4925 
4926  INT n_coarse = row;
4927  INT n_fine = A->ROW;
4928  INT coarse_mul_nthreads = n_coarse * nthreads;
4929  INT fine_mul_nthreads = n_fine * nthreads;
4930  INT coarse_add_nthreads = n_coarse + nthreads;
4931  INT minus_one_length = coarse_mul_nthreads + fine_mul_nthreads;
4932  INT total_calloc = minus_one_length + coarse_add_nthreads + nthreads;
4933 
4934  Ps_marker = (INT *)fasp_mem_calloc(total_calloc, sizeof(INT));
4935  As_marker = Ps_marker + coarse_mul_nthreads;
4936 
4937  /*------------------------------------------------------*
4938  * First Pass: Determine size of B and set up B_i *
4939  *------------------------------------------------------*/
4940  iac = (INT *)fasp_mem_calloc(n_coarse+1, sizeof(INT));
4941 
4942  fasp_iarray_set(minus_one_length, Ps_marker, -1);
4943 
4944  REAL *tmp=(REAL *)fasp_mem_calloc(2*nthreads*nb2, sizeof(REAL));
4945 
4946 #ifdef _OPENMP
4947  INT * RAP_temp = As_marker + fine_mul_nthreads;
4948  INT * part_end = RAP_temp + coarse_add_nthreads;
4949 
4950  if (n_coarse > OPENMP_HOLDS) {
4951 #pragma omp parallel for private(myid, mybegin, myend, Ctemp, P_marker, A_marker, \
4952 counter, i, jj_row_begining, jj1, i1, jj2, i2, jj3, i3)
4953  for (myid = 0; myid < nthreads; myid++) {
4954  FASP_GET_START_END(myid, nthreads, n_coarse, &mybegin, &myend);
4955  P_marker = Ps_marker + myid * n_coarse;
4956  A_marker = As_marker + myid * n_fine;
4957  counter = 0;
4958  for (i = mybegin; i < myend; ++i) {
4959  P_marker[i] = counter;
4960  jj_row_begining = counter;
4961  counter ++;
4962  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
4963  i1 = jr[jj1];
4964  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
4965  i2 = ja[jj2];
4966  if (A_marker[i2] != i) {
4967  A_marker[i2] = i;
4968  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
4969  i3 = jp[jj3];
4970  if (P_marker[i3] < jj_row_begining) {
4971  P_marker[i3] = counter;
4972  counter ++;
4973  }
4974  }
4975  }
4976  }
4977  }
4978  RAP_temp[i+myid] = jj_row_begining;
4979  }
4980  RAP_temp[myend+myid] = counter;
4981  part_end[myid] = myend + myid + 1;
4982  }
4983  fasp_iarray_cp(part_end[0], RAP_temp, iac);
4984  counter = part_end[0];
4985  Ctemp = 0;
4986  for (i1 = 1; i1 < nthreads; i1 ++) {
4987  Ctemp += RAP_temp[part_end[i1-1]-1];
4988  for (jj1 = part_end[i1-1]+1; jj1 < part_end[i1]; jj1 ++) {
4989  iac[counter] = RAP_temp[jj1] + Ctemp;
4990  counter ++;
4991  }
4992  }
4993  }
4994  else {
4995 #endif
4996  counter = 0;
4997  for (i = 0; i < row; ++ i) {
4998  Ps_marker[i] = counter;
4999  jj_row_begining = counter;
5000  counter ++;
5001 
5002  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5003  i1 = jr[jj1];
5004  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5005  i2 = ja[jj2];
5006  if (As_marker[i2] != i) {
5007  As_marker[i2] = i;
5008  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5009  i3 = jp[jj3];
5010  if (Ps_marker[i3] < jj_row_begining) {
5011  Ps_marker[i3] = counter;
5012  counter ++;
5013  }
5014  }
5015  }
5016  }
5017  }
5018  iac[i] = jj_row_begining;
5019  }
5020 #ifdef _OPENMP
5021  }
5022 #endif
5023 
5024  iac[row] = counter;
5025 
5026  jac=(INT*)fasp_mem_calloc(iac[row], sizeof(INT));
5027 
5028  acj=(REAL*)fasp_mem_calloc(iac[row]*nb2, sizeof(REAL));
5029 
5030  fasp_iarray_set(minus_one_length, Ps_marker, -1);
5031 
5032  /*------------------------------------------------------*
5033  * Second Pass: compute entries of B=R*A*P *
5034  *------------------------------------------------------*/
5035 #ifdef _OPENMP
5036  if (n_coarse > OPENMP_HOLDS) {
5037 #pragma omp parallel for private(myid, mybegin, myend, Ctemp, P_marker, A_marker, \
5038 counter, i, jj_row_begining, jj1, i1, jj2, i2, \
5039 jj3, i3, smat_tmp)
5040  for (myid = 0; myid < nthreads; myid++) {
5041  FASP_GET_START_END(myid, nthreads, n_coarse, &mybegin, &myend);
5042  P_marker = Ps_marker + myid * n_coarse;
5043  A_marker = As_marker + myid * n_fine;
5044  smat_tmp = tmp + myid*2*nb2;
5045  counter = iac[mybegin];
5046  for (i = mybegin; i < myend; ++i) {
5047  P_marker[i] = counter;
5048  jj_row_begining = counter;
5049  jac[counter] = i;
5050  fasp_array_set(nb2, &acj[counter*nb2], 0x0);
5051  counter ++;
5052 
5053  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5054  i1 = jr[jj1];
5055  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5056  fasp_blas_smat_mul(&rj[jj1*nb2],&aj[jj2*nb2], smat_tmp, nb);
5057  i2 = ja[jj2];
5058  if (A_marker[i2] != i) {
5059  A_marker[i2] = i;
5060  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5061  i3 = jp[jj3];
5062  fasp_blas_smat_mul(smat_tmp, &pj[jj3*nb2], smat_tmp+nb2, nb);
5063  if (P_marker[i3] < jj_row_begining) {
5064  P_marker[i3] = counter;
5065  fasp_array_cp(nb2, smat_tmp+nb2, &acj[counter*nb2]);
5066  jac[counter] = i3;
5067  counter ++;
5068  }
5069  else {
5070  fasp_blas_array_axpy(nb2, 1.0, smat_tmp+nb2,
5071  &acj[P_marker[i3]*nb2]);
5072  }
5073  }
5074  }
5075  else {
5076  for (jj3 = ip[i2]; jj3 < ip[i2+1]; jj3 ++) {
5077  i3 = jp[jj3];
5078  fasp_blas_smat_mul(smat_tmp, &pj[jj3*nb2], smat_tmp+nb2, nb);
5079  fasp_blas_array_axpy(nb2, 1.0, smat_tmp+nb2,
5080  &acj[P_marker[i3]*nb2]);
5081  }
5082  }
5083  }
5084  }
5085  }
5086  }
5087  }
5088  else {
5089 #endif
5090  counter = 0;
5091  for (i = 0; i < row; ++i) {
5092  Ps_marker[i] = counter;
5093  jj_row_begining = counter;
5094  jac[counter] = i;
5095  fasp_array_set(nb2, &acj[counter*nb2], 0x0);
5096  counter ++;
5097 
5098  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5099  i1 = jr[jj1];
5100  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5101  fasp_blas_smat_mul(&rj[jj1*nb2],&aj[jj2*nb2], tmp, nb);
5102  i2 = ja[jj2];
5103  if (As_marker[i2] != i) {
5104  As_marker[i2] = i;
5105  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5106  i3 = jp[jj3];
5107  fasp_blas_smat_mul(tmp, &pj[jj3*nb2], tmp+nb2, nb);
5108  if (Ps_marker[i3] < jj_row_begining) {
5109  Ps_marker[i3] = counter;
5110  fasp_array_cp(nb2, tmp+nb2, &acj[counter*nb2]);
5111  jac[counter] = i3;
5112  counter ++;
5113  }
5114  else {
5115  fasp_blas_array_axpy(nb2, 1.0, tmp+nb2,
5116  &acj[Ps_marker[i3]*nb2]);
5117  }
5118  }
5119  }
5120  else {
5121  for (jj3 = ip[i2]; jj3 < ip[i2+1]; jj3 ++) {
5122  i3 = jp[jj3];
5123  fasp_blas_smat_mul(tmp, &pj[jj3*nb2], tmp+nb2, nb);
5124  fasp_blas_array_axpy(nb2, 1.0, tmp+nb2, &acj[Ps_marker[i3]*nb2]);
5125  }
5126  }
5127  }
5128  }
5129  }
5130 #ifdef _OPENMP
5131  }
5132 #endif
5133  // setup coarse matrix B
5134  B->ROW=row; B->COL=col;
5135  B->IA=iac; B->JA=jac; B->val=acj;
5136  B->NNZ=B->IA[B->ROW]-B->IA[0];
5137  B->nb=A->nb;
5139  if(Ps_marker) fasp_mem_free(Ps_marker);
5140  if(tmp) fasp_mem_free(tmp);
5141 }
5142 
5161  dBSRmat *A,
5162  dBSRmat *P,
5163  dBSRmat *B)
5164 {
5165  const INT row=R->ROW, col=P->COL, nb2=A->nb*A->nb;
5166  //INT nB=A->NNZ;
5167 
5168  REAL *aj=A->val, *acj;
5169  INT *ir=R->IA, *ia=A->IA, *ip=P->IA, *iac;
5170  INT *jr=R->JA, *ja=A->JA, *jp=P->JA, *jac;
5171 
5172  INT *Ps_marker = NULL;
5173  INT *As_marker = NULL;
5174 
5175 #ifdef _OPENMP
5176  INT *P_marker = NULL;
5177  INT *A_marker = NULL;
5178  REAL *smat_tmp = NULL;
5179 #endif
5180 
5181  INT i, i1, i2, i3, jj1, jj2, jj3;
5182  INT counter, jj_row_begining;
5183 
5184  INT nthreads = 1;
5185 
5186 #ifdef _OPENMP
5187  INT myid, mybegin, myend, Ctemp;
5188  nthreads = FASP_GET_NUM_THREADS();
5189 #endif
5190 
5191  INT n_coarse = row;
5192  INT n_fine = A->ROW;
5193  INT coarse_mul_nthreads = n_coarse * nthreads;
5194  INT fine_mul_nthreads = n_fine * nthreads;
5195  INT coarse_add_nthreads = n_coarse + nthreads;
5196  INT minus_one_length = coarse_mul_nthreads + fine_mul_nthreads;
5197  INT total_calloc = minus_one_length + coarse_add_nthreads + nthreads;
5198 
5199  Ps_marker = (INT *)fasp_mem_calloc(total_calloc, sizeof(INT));
5200  As_marker = Ps_marker + coarse_mul_nthreads;
5201 
5202  /*------------------------------------------------------*
5203  * First Pass: Determine size of B and set up B_i *
5204  *------------------------------------------------------*/
5205  iac = (INT *)fasp_mem_calloc(n_coarse+1, sizeof(INT));
5206 
5207  fasp_iarray_set(minus_one_length, Ps_marker, -1);
5208 
5209  REAL *tmp=(REAL *)fasp_mem_calloc(2*nthreads*nb2, sizeof(REAL));
5210 
5211 #ifdef _OPENMP
5212  INT * RAP_temp = As_marker + fine_mul_nthreads;
5213  INT * part_end = RAP_temp + coarse_add_nthreads;
5214 
5215  if (n_coarse > OPENMP_HOLDS) {
5216 #pragma omp parallel for private(myid, mybegin, myend, Ctemp, P_marker, A_marker, \
5217 counter, i, jj_row_begining, jj1, i1, jj2, i2, jj3, i3)
5218  for (myid = 0; myid < nthreads; myid++) {
5219  FASP_GET_START_END(myid, nthreads, n_coarse, &mybegin, &myend);
5220  P_marker = Ps_marker + myid * n_coarse;
5221  A_marker = As_marker + myid * n_fine;
5222  counter = 0;
5223  for (i = mybegin; i < myend; ++i) {
5224  P_marker[i] = counter;
5225  jj_row_begining = counter;
5226  counter ++;
5227  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5228  i1 = jr[jj1];
5229  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5230  i2 = ja[jj2];
5231  if (A_marker[i2] != i) {
5232  A_marker[i2] = i;
5233  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5234  i3 = jp[jj3];
5235  if (P_marker[i3] < jj_row_begining) {
5236  P_marker[i3] = counter;
5237  counter ++;
5238  }
5239  }
5240  }
5241  }
5242  }
5243  RAP_temp[i+myid] = jj_row_begining;
5244  }
5245  RAP_temp[myend+myid] = counter;
5246  part_end[myid] = myend + myid + 1;
5247  }
5248  fasp_iarray_cp(part_end[0], RAP_temp, iac);
5249  counter = part_end[0];
5250  Ctemp = 0;
5251  for (i1 = 1; i1 < nthreads; i1 ++) {
5252  Ctemp += RAP_temp[part_end[i1-1]-1];
5253  for (jj1 = part_end[i1-1]+1; jj1 < part_end[i1]; jj1 ++) {
5254  iac[counter] = RAP_temp[jj1] + Ctemp;
5255  counter ++;
5256  }
5257  }
5258  }
5259  else {
5260 #endif
5261  counter = 0;
5262  for (i = 0; i < row; ++ i) {
5263  Ps_marker[i] = counter;
5264  jj_row_begining = counter;
5265  counter ++;
5266 
5267  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5268  i1 = jr[jj1];
5269  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5270  i2 = ja[jj2];
5271  if (As_marker[i2] != i) {
5272  As_marker[i2] = i;
5273  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5274  i3 = jp[jj3];
5275  if (Ps_marker[i3] < jj_row_begining) {
5276  Ps_marker[i3] = counter;
5277  counter ++;
5278  }
5279  }
5280  }
5281  }
5282  }
5283  iac[i] = jj_row_begining;
5284  }
5285 #ifdef _OPENMP
5286  }
5287 #endif
5288 
5289  iac[row] = counter;
5290 
5291  jac=(INT*)fasp_mem_calloc(iac[row], sizeof(INT));
5292 
5293  acj=(REAL*)fasp_mem_calloc(iac[row]*nb2, sizeof(REAL));
5294 
5295  fasp_iarray_set(minus_one_length, Ps_marker, -1);
5296 
5297  /*------------------------------------------------------*
5298  * Second Pass: compute entries of B=R*A*P *
5299  *------------------------------------------------------*/
5300 #ifdef _OPENMP
5301  if (n_coarse > OPENMP_HOLDS) {
5302 #pragma omp parallel for private(myid, mybegin, myend, Ctemp, P_marker, A_marker, \
5303 counter, i, jj_row_begining, jj1, i1, jj2, i2, \
5304 jj3, i3, smat_tmp)
5305  for (myid = 0; myid < nthreads; myid++) {
5306  FASP_GET_START_END(myid, nthreads, n_coarse, &mybegin, &myend);
5307  P_marker = Ps_marker + myid * n_coarse;
5308  A_marker = As_marker + myid * n_fine;
5309  smat_tmp = tmp + myid*2*nb2;
5310  counter = iac[mybegin];
5311  for (i = mybegin; i < myend; ++i) {
5312  P_marker[i] = counter;
5313  jj_row_begining = counter;
5314  jac[counter] = i;
5315  fasp_array_set(nb2, &acj[counter*nb2], 0x0);
5316  counter ++;
5317 
5318  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5319  i1 = jr[jj1];
5320  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5321 
5322  i2 = ja[jj2];
5323  if (A_marker[i2] != i) {
5324  A_marker[i2] = i;
5325  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5326  i3 = jp[jj3];
5327 
5328  if (P_marker[i3] < jj_row_begining) {
5329  P_marker[i3] = counter;
5330  fasp_array_cp(nb2, &aj[jj2*nb2], &acj[counter*nb2]);
5331  jac[counter] = i3;
5332  counter ++;
5333  }
5334  else {
5335  fasp_blas_array_axpy(nb2, 1.0, &aj[jj2*nb2],
5336  &acj[P_marker[i3]*nb2]);
5337  }
5338  }
5339  }
5340  else {
5341  for (jj3 = ip[i2]; jj3 < ip[i2+1]; jj3 ++) {
5342  i3 = jp[jj3];
5343  fasp_blas_array_axpy(nb2, 1.0, &aj[jj2*nb2],
5344  &acj[P_marker[i3]*nb2]);
5345  }
5346  }
5347  }
5348  }
5349  }
5350  }
5351  }
5352  else {
5353 #endif
5354  counter = 0;
5355  for (i = 0; i < row; ++i) {
5356  Ps_marker[i] = counter;
5357  jj_row_begining = counter;
5358  jac[counter] = i;
5359  fasp_array_set(nb2, &acj[counter*nb2], 0x0);
5360  counter ++;
5361 
5362  for (jj1 = ir[i]; jj1 < ir[i+1]; ++jj1) {
5363  i1 = jr[jj1];
5364  for (jj2 = ia[i1]; jj2 < ia[i1+1]; ++jj2) {
5365 
5366  i2 = ja[jj2];
5367  if (As_marker[i2] != i) {
5368  As_marker[i2] = i;
5369  for (jj3 = ip[i2]; jj3 < ip[i2+1]; ++jj3) {
5370  i3 = jp[jj3];
5371  if (Ps_marker[i3] < jj_row_begining) {
5372  Ps_marker[i3] = counter;
5373  fasp_array_cp(nb2, &aj[jj2*nb2], &acj[counter*nb2]);
5374  jac[counter] = i3;
5375  counter ++;
5376  }
5377  else {
5378  fasp_blas_array_axpy(nb2, 1.0, &aj[jj2*nb2],
5379  &acj[Ps_marker[i3]*nb2]);
5380  }
5381  }
5382  }
5383  else {
5384  for (jj3 = ip[i2]; jj3 < ip[i2+1]; jj3 ++) {
5385  i3 = jp[jj3];
5386  fasp_blas_array_axpy(nb2, 1.0, &aj[jj2*nb2],
5387  &acj[Ps_marker[i3]*nb2]);
5388  }
5389  }
5390  }
5391  }
5392  }
5393 #ifdef _OPENMP
5394  }
5395 #endif
5396 
5397  // setup coarse matrix B
5398  B->ROW=row; B->COL=col;
5399  B->IA=iac; B->JA=jac; B->val=acj;
5400  B->NNZ=B->IA[B->ROW]-B->IA[0];
5401  B->nb=A->nb;
5403  if(Ps_marker) fasp_mem_free(Ps_marker);
5404  if(tmp) fasp_mem_free(tmp);
5405 }
5406 
5407 /*---------------------------------*/
5408 /*-- End of File --*/
5409 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
void fasp_blas_smat_ypAx(REAL *A, REAL *x, REAL *y, const INT n)
Compute y := y + Ax, where 'A' is a n*n dense matrix.
Definition: blas_smat.c:976
#define REAL
Definition: fasp.h:67
void fasp_blas_dbsr_mxm(dBSRmat *A, dBSRmat *B, dBSRmat *C)
Sparse matrix multiplication C=A*B.
Definition: blas_bsr.c:4591
INT nb
dimension of each sub-block
Definition: fasp_block.h:56
void fasp_blas_dbsr_aAxpy_agg(const REAL alpha, dBSRmat *A, REAL *x, REAL *y)
Compute y := alpha*A*x + y where each small block matrix is an identity matrix.
Definition: blas_bsr.c:610
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
void fasp_blas_smat_ypAx_nc7(REAL *A, REAL *x, REAL *y)
Compute y := y + Ax, where 'A' is a 7*7 dense matrix.
Definition: blas_smat.c:941
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void fasp_blas_dbsr_rap_agg(dBSRmat *R, dBSRmat *A, dBSRmat *P, dBSRmat *B)
dBSRmat sparse matrix multiplication B=R*A*P, where small block matrices in P and R are identity matr...
Definition: blas_bsr.c:5160
void fasp_blas_smat_mul(REAL *a, REAL *b, REAL *c, const INT n)
Compute the matrix product of two small full matrices a and b, stored in c.
Definition: blas_smat.c:448
#define INT
Definition: fasp.h:64
void fasp_blas_smat_ypAx_nc3(REAL *A, REAL *x, REAL *y)
Compute y := y + Ax, where 'A' is a 3*3 dense matrix.
Definition: blas_smat.c:883
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
REAL * val
Definition: fasp_block.h:67
INT NNZ
number of nonzero sub-blocks in matrix A, NNZ
Definition: fasp_block.h:53
void fasp_iarray_cp(const INT n, INT *x, INT *y)
Copy an array to the other y=x.
Definition: array.c:185
void fasp_blas_dbsr_aAxpby(const REAL alpha, dBSRmat *A, REAL *x, const REAL beta, REAL *y)
Compute y := alpha*A*x + beta*y.
Definition: blas_bsr.c:59
#define OPENMP_HOLDS
Definition: fasp_const.h:248
INT storage_manner
storage manner for each sub-block
Definition: fasp_block.h:59
Main header file for FASP.
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
void fasp_blas_smat_ypAx_nc2(REAL *A, REAL *x, REAL *y)
Compute y := y + Ax, where 'A' is a 2*2 dense matrix.
Definition: blas_smat.c:857
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
INT COL
number of cols of sub-blocks in matrix A, N
Definition: fasp_block.h:50
void fasp_blas_dbsr_mxv(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x.
Definition: blas_bsr.c:895
void fasp_blas_smat_ypAx_nc5(REAL *A, REAL *x, REAL *y)
Compute y := y + Ax, where 'A' is a 5*5 dense matrix.
Definition: blas_smat.c:910
INT count
void fasp_blas_dbsr_rap(dBSRmat *R, dBSRmat *A, dBSRmat *P, dBSRmat *B)
dBSRmat sparse matrix multiplication B=R*A*P
Definition: blas_bsr.c:4895
void fasp_blas_array_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: blas_array.c:35
void fasp_blas_dbsr_rap1(dBSRmat *R, dBSRmat *A, dBSRmat *P, dBSRmat *B)
dBSRmat sparse matrix multiplication B=R*A*P
Definition: blas_bsr.c:4711
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:47
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
void fasp_blas_dbsr_mxv_agg(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x, where each small block matrices of A is an identity matrix.
Definition: blas_bsr.c:2641
INT * JA
Definition: fasp_block.h:74
#define FALSE
Definition: fasp_const.h:68
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: array.c:107
INT * IA
integer array of row pointers, the size is ROW+1
Definition: fasp_block.h:70
void fasp_blas_dbsr_axm(dBSRmat *A, const REAL alpha)
Multiply a sparse matrix A in BSR format by a scalar alpha.
Definition: blas_bsr.c:30
void fasp_blas_dbsr_aAxpy(const REAL alpha, dBSRmat *A, REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: blas_bsr.c:337