← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

Viewing changes to src/mat/impls/maij/maij.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mat/impls/maij/maij.c
17
17
*/
18
18
 
19
19
#include "src/mat/impls/maij/maij.h"
 
20
#include "vecimpl.h"
20
21
 
21
22
#undef __FUNCT__  
22
23
#define __FUNCT__ "MatMAIJGetAIJ" 
97
98
  PetscFunctionReturn(0);
98
99
}
99
100
 
 
101
/*MC
 
102
  MATMAIJ - MATMAIJ = "maij" - A matrix type to be used for restriction and interpolation operations for 
 
103
  multicomponent problems, interpolating or restricting each component the same way independently.
 
104
  The matrix type is based on MATSEQAIJ for sequential matrices, and MATMPIAIJ for distributed matrices.
 
105
 
 
106
  Operations provided:
 
107
. MatMult
 
108
. MatMultTranspose
 
109
. MatMultAdd
 
110
. MatMultTransposeAdd
 
111
 
 
112
  Level: advanced
 
113
 
 
114
.seealso: MatCreateSeqDense
 
115
M*/
 
116
 
100
117
EXTERN_C_BEGIN
101
118
#undef __FUNCT__  
102
119
#define __FUNCT__ "MatCreate_MAIJ" 
130
147
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
131
148
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
132
149
  PetscScalar   *x,*y,*v,sum1, sum2;
133
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
150
  int           ierr,m = b->AIJ->m,*idx,*ii;
134
151
  int           n,i,jrow,j;
135
152
 
136
153
  PetscFunctionBegin;
137
154
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
138
155
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
139
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
140
156
  idx  = a->j;
141
157
  v    = a->a;
142
158
  ii   = a->i;
143
159
 
144
 
  v    += shift; /* shift for Fortran start by 1 indexing */
145
 
  idx  += shift;
146
160
  for (i=0; i<m; i++) {
147
161
    jrow = ii[i];
148
162
    n    = ii[i+1] - jrow;
170
184
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
171
185
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
172
186
  PetscScalar   *x,*y,*v,alpha1,alpha2,zero = 0.0;
173
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
187
  int           ierr,m = b->AIJ->m,n,i,*idx;
174
188
 
175
189
  PetscFunctionBegin; 
176
190
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
177
191
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
178
192
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
179
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
193
 
180
194
  for (i=0; i<m; i++) {
181
 
    idx    = a->j + a->i[i] + shift;
182
 
    v      = a->a + a->i[i] + shift;
 
195
    idx    = a->j + a->i[i] ;
 
196
    v      = a->a + a->i[i] ;
183
197
    n      = a->i[i+1] - a->i[i];
184
198
    alpha1 = x[2*i];
185
199
    alpha2 = x[2*i+1];
198
212
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
199
213
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
200
214
  PetscScalar   *x,*y,*v,sum1, sum2;
201
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
215
  int           ierr,m = b->AIJ->m,*idx,*ii;
202
216
  int           n,i,jrow,j;
203
217
 
204
218
  PetscFunctionBegin;
205
219
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
206
220
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
207
221
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
208
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
209
222
  idx  = a->j;
210
223
  v    = a->a;
211
224
  ii   = a->i;
212
225
 
213
 
  v    += shift; /* shift for Fortran start by 1 indexing */
214
 
  idx  += shift;
215
226
  for (i=0; i<m; i++) {
216
227
    jrow = ii[i];
217
228
    n    = ii[i+1] - jrow;
238
249
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
239
250
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
240
251
  PetscScalar   *x,*y,*v,alpha1,alpha2;
241
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
252
  int           ierr,m = b->AIJ->m,n,i,*idx;
242
253
 
243
254
  PetscFunctionBegin; 
244
255
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
245
256
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
246
257
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
247
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
258
 
248
259
  for (i=0; i<m; i++) {
249
 
    idx   = a->j + a->i[i] + shift;
250
 
    v     = a->a + a->i[i] + shift;
 
260
    idx   = a->j + a->i[i] ;
 
261
    v     = a->a + a->i[i] ;
251
262
    n     = a->i[i+1] - a->i[i];
252
263
    alpha1 = x[2*i];
253
264
    alpha2 = x[2*i+1];
266
277
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
267
278
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
268
279
  PetscScalar   *x,*y,*v,sum1, sum2, sum3;
269
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
280
  int           ierr,m = b->AIJ->m,*idx,*ii;
270
281
  int           n,i,jrow,j;
271
282
 
272
283
  PetscFunctionBegin;
273
284
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
274
285
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
275
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
276
286
  idx  = a->j;
277
287
  v    = a->a;
278
288
  ii   = a->i;
279
289
 
280
 
  v    += shift; /* shift for Fortran start by 1 indexing */
281
 
  idx  += shift;
282
290
  for (i=0; i<m; i++) {
283
291
    jrow = ii[i];
284
292
    n    = ii[i+1] - jrow;
309
317
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
310
318
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
311
319
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,zero = 0.0;
312
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
320
  int           ierr,m = b->AIJ->m,n,i,*idx;
313
321
 
314
322
  PetscFunctionBegin; 
315
323
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
316
324
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
317
325
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
318
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
326
  
319
327
  for (i=0; i<m; i++) {
320
 
    idx    = a->j + a->i[i] + shift;
321
 
    v      = a->a + a->i[i] + shift;
 
328
    idx    = a->j + a->i[i];
 
329
    v      = a->a + a->i[i];
322
330
    n      = a->i[i+1] - a->i[i];
323
331
    alpha1 = x[3*i];
324
332
    alpha2 = x[3*i+1];
343
351
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
344
352
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
345
353
  PetscScalar   *x,*y,*v,sum1, sum2, sum3;
346
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
354
  int           ierr,m = b->AIJ->m,*idx,*ii;
347
355
  int           n,i,jrow,j;
348
356
 
349
357
  PetscFunctionBegin;
350
358
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
351
359
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
352
360
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
353
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
354
361
  idx  = a->j;
355
362
  v    = a->a;
356
363
  ii   = a->i;
357
364
 
358
 
  v    += shift; /* shift for Fortran start by 1 indexing */
359
 
  idx  += shift;
360
365
  for (i=0; i<m; i++) {
361
366
    jrow = ii[i];
362
367
    n    = ii[i+1] - jrow;
386
391
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
387
392
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
388
393
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3;
389
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
394
  int           ierr,m = b->AIJ->m,n,i,*idx;
390
395
 
391
396
  PetscFunctionBegin; 
392
397
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
393
398
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
394
399
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
395
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
396
400
  for (i=0; i<m; i++) {
397
 
    idx    = a->j + a->i[i] + shift;
398
 
    v      = a->a + a->i[i] + shift;
 
401
    idx    = a->j + a->i[i] ;
 
402
    v      = a->a + a->i[i] ;
399
403
    n      = a->i[i+1] - a->i[i];
400
404
    alpha1 = x[3*i];
401
405
    alpha2 = x[3*i+1];
421
425
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
422
426
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
423
427
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4;
424
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
428
  int           ierr,m = b->AIJ->m,*idx,*ii;
425
429
  int           n,i,jrow,j;
426
430
 
427
431
  PetscFunctionBegin;
428
432
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
429
433
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
430
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
431
434
  idx  = a->j;
432
435
  v    = a->a;
433
436
  ii   = a->i;
434
437
 
435
 
  v    += shift; /* shift for Fortran start by 1 indexing */
436
 
  idx  += shift;
437
438
  for (i=0; i<m; i++) {
438
439
    jrow = ii[i];
439
440
    n    = ii[i+1] - jrow;
467
468
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
468
469
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
469
470
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,zero = 0.0;
470
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
471
  int           ierr,m = b->AIJ->m,n,i,*idx;
471
472
 
472
473
  PetscFunctionBegin; 
473
474
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
474
475
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
475
476
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
476
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
477
477
  for (i=0; i<m; i++) {
478
 
    idx    = a->j + a->i[i] + shift;
479
 
    v      = a->a + a->i[i] + shift;
 
478
    idx    = a->j + a->i[i] ;
 
479
    v      = a->a + a->i[i] ;
480
480
    n      = a->i[i+1] - a->i[i];
481
481
    alpha1 = x[4*i];
482
482
    alpha2 = x[4*i+1];
503
503
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
504
504
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
505
505
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4;
506
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
506
  int           ierr,m = b->AIJ->m,*idx,*ii;
507
507
  int           n,i,jrow,j;
508
508
 
509
509
  PetscFunctionBegin;
510
510
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
511
511
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
512
512
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
513
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
514
513
  idx  = a->j;
515
514
  v    = a->a;
516
515
  ii   = a->i;
517
516
 
518
 
  v    += shift; /* shift for Fortran start by 1 indexing */
519
 
  idx  += shift;
520
517
  for (i=0; i<m; i++) {
521
518
    jrow = ii[i];
522
519
    n    = ii[i+1] - jrow;
549
546
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
550
547
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
551
548
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4;
552
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
549
  int           ierr,m = b->AIJ->m,n,i,*idx;
553
550
 
554
551
  PetscFunctionBegin; 
555
552
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
556
553
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
557
554
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
558
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
555
  
559
556
  for (i=0; i<m; i++) {
560
 
    idx    = a->j + a->i[i] + shift;
561
 
    v      = a->a + a->i[i] + shift;
 
557
    idx    = a->j + a->i[i] ;
 
558
    v      = a->a + a->i[i] ;
562
559
    n      = a->i[i+1] - a->i[i];
563
560
    alpha1 = x[4*i];
564
561
    alpha2 = x[4*i+1];
586
583
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
587
584
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
588
585
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5;
589
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
586
  int           ierr,m = b->AIJ->m,*idx,*ii;
590
587
  int           n,i,jrow,j;
591
588
 
592
589
  PetscFunctionBegin;
593
590
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
594
591
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
595
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
596
592
  idx  = a->j;
597
593
  v    = a->a;
598
594
  ii   = a->i;
599
595
 
600
 
  v    += shift; /* shift for Fortran start by 1 indexing */
601
 
  idx  += shift;
602
596
  for (i=0; i<m; i++) {
603
597
    jrow = ii[i];
604
598
    n    = ii[i+1] - jrow;
635
629
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
636
630
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
637
631
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,zero = 0.0;
638
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
632
  int           ierr,m = b->AIJ->m,n,i,*idx;
639
633
 
640
634
  PetscFunctionBegin; 
641
635
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
642
636
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
643
637
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
644
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
638
  
645
639
  for (i=0; i<m; i++) {
646
 
    idx    = a->j + a->i[i] + shift;
647
 
    v      = a->a + a->i[i] + shift;
 
640
    idx    = a->j + a->i[i] ;
 
641
    v      = a->a + a->i[i] ;
648
642
    n      = a->i[i+1] - a->i[i];
649
643
    alpha1 = x[5*i];
650
644
    alpha2 = x[5*i+1];
673
667
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
674
668
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
675
669
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5;
676
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
670
  int           ierr,m = b->AIJ->m,*idx,*ii;
677
671
  int           n,i,jrow,j;
678
672
 
679
673
  PetscFunctionBegin;
680
674
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
681
675
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
682
676
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
683
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
684
677
  idx  = a->j;
685
678
  v    = a->a;
686
679
  ii   = a->i;
687
680
 
688
 
  v    += shift; /* shift for Fortran start by 1 indexing */
689
 
  idx  += shift;
690
681
  for (i=0; i<m; i++) {
691
682
    jrow = ii[i];
692
683
    n    = ii[i+1] - jrow;
723
714
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
724
715
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
725
716
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5;
726
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
717
  int           ierr,m = b->AIJ->m,n,i,*idx;
727
718
 
728
719
  PetscFunctionBegin; 
729
720
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
730
721
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
731
722
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
732
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
723
 
733
724
  for (i=0; i<m; i++) {
734
 
    idx    = a->j + a->i[i] + shift;
735
 
    v      = a->a + a->i[i] + shift;
 
725
    idx    = a->j + a->i[i] ;
 
726
    v      = a->a + a->i[i] ;
736
727
    n      = a->i[i+1] - a->i[i];
737
728
    alpha1 = x[5*i];
738
729
    alpha2 = x[5*i+1];
762
753
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
763
754
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
764
755
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6;
765
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
756
  int           ierr,m = b->AIJ->m,*idx,*ii;
766
757
  int           n,i,jrow,j;
767
758
 
768
759
  PetscFunctionBegin;
769
760
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
770
761
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
771
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
772
762
  idx  = a->j;
773
763
  v    = a->a;
774
764
  ii   = a->i;
775
765
 
776
 
  v    += shift; /* shift for Fortran start by 1 indexing */
777
 
  idx  += shift;
778
766
  for (i=0; i<m; i++) {
779
767
    jrow = ii[i];
780
768
    n    = ii[i+1] - jrow;
814
802
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
815
803
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
816
804
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,zero = 0.0;
817
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
805
  int           ierr,m = b->AIJ->m,n,i,*idx;
818
806
 
819
807
  PetscFunctionBegin; 
820
808
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
821
809
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
822
810
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
823
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
811
 
824
812
  for (i=0; i<m; i++) {
825
 
    idx    = a->j + a->i[i] + shift;
826
 
    v      = a->a + a->i[i] + shift;
 
813
    idx    = a->j + a->i[i] ;
 
814
    v      = a->a + a->i[i] ;
827
815
    n      = a->i[i+1] - a->i[i];
828
816
    alpha1 = x[6*i];
829
817
    alpha2 = x[6*i+1];
854
842
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
855
843
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
856
844
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6;
857
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
845
  int           ierr,m = b->AIJ->m,*idx,*ii;
858
846
  int           n,i,jrow,j;
859
847
 
860
848
  PetscFunctionBegin;
861
849
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
862
850
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
863
851
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
864
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
865
852
  idx  = a->j;
866
853
  v    = a->a;
867
854
  ii   = a->i;
868
855
 
869
 
  v    += shift; /* shift for Fortran start by 1 indexing */
870
 
  idx  += shift;
871
856
  for (i=0; i<m; i++) {
872
857
    jrow = ii[i];
873
858
    n    = ii[i+1] - jrow;
907
892
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
908
893
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
909
894
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6;
910
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
895
  int           ierr,m = b->AIJ->m,n,i,*idx;
911
896
 
912
897
  PetscFunctionBegin; 
913
898
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
914
899
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
915
900
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
916
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
901
 
917
902
  for (i=0; i<m; i++) {
918
 
    idx    = a->j + a->i[i] + shift;
919
 
    v      = a->a + a->i[i] + shift;
 
903
    idx    = a->j + a->i[i] ;
 
904
    v      = a->a + a->i[i] ;
920
905
    n      = a->i[i+1] - a->i[i];
921
906
    alpha1 = x[6*i];
922
907
    alpha2 = x[6*i+1];
948
933
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
949
934
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
950
935
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8;
951
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
936
  int           ierr,m = b->AIJ->m,*idx,*ii;
952
937
  int           n,i,jrow,j;
953
938
 
954
939
  PetscFunctionBegin;
955
940
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
956
941
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
957
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
958
942
  idx  = a->j;
959
943
  v    = a->a;
960
944
  ii   = a->i;
961
945
 
962
 
  v    += shift; /* shift for Fortran start by 1 indexing */
963
 
  idx  += shift;
964
946
  for (i=0; i<m; i++) {
965
947
    jrow = ii[i];
966
948
    n    = ii[i+1] - jrow;
1006
988
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
1007
989
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
1008
990
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8,zero = 0.0;
1009
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
991
  int           ierr,m = b->AIJ->m,n,i,*idx;
1010
992
 
1011
993
  PetscFunctionBegin; 
1012
994
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
1013
995
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1014
996
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
1015
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
 
997
 
1016
998
  for (i=0; i<m; i++) {
1017
 
    idx    = a->j + a->i[i] + shift;
1018
 
    v      = a->a + a->i[i] + shift;
 
999
    idx    = a->j + a->i[i] ;
 
1000
    v      = a->a + a->i[i] ;
1019
1001
    n      = a->i[i+1] - a->i[i];
1020
1002
    alpha1 = x[8*i];
1021
1003
    alpha2 = x[8*i+1];
1050
1032
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
1051
1033
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
1052
1034
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8;
1053
 
  int           ierr,m = b->AIJ->m,*idx,shift = a->indexshift,*ii;
 
1035
  int           ierr,m = b->AIJ->m,*idx,*ii;
1054
1036
  int           n,i,jrow,j;
1055
1037
 
1056
1038
  PetscFunctionBegin;
1057
1039
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
1058
1040
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1059
1041
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
1060
 
  x    = x + shift;    /* shift for Fortran start by 1 indexing */
1061
1042
  idx  = a->j;
1062
1043
  v    = a->a;
1063
1044
  ii   = a->i;
1064
1045
 
1065
 
  v    += shift; /* shift for Fortran start by 1 indexing */
1066
 
  idx  += shift;
1067
1046
  for (i=0; i<m; i++) {
1068
1047
    jrow = ii[i];
1069
1048
    n    = ii[i+1] - jrow;
1109
1088
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
1110
1089
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
1111
1090
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8;
1112
 
  int           ierr,m = b->AIJ->m,n,i,*idx,shift = a->indexshift;
 
1091
  int           ierr,m = b->AIJ->m,n,i,*idx;
1113
1092
 
1114
1093
  PetscFunctionBegin; 
1115
1094
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
1116
1095
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1117
1096
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
1118
 
  y = y + shift; /* shift for Fortran start by 1 indexing */
1119
1097
  for (i=0; i<m; i++) {
1120
 
    idx    = a->j + a->i[i] + shift;
1121
 
    v      = a->a + a->i[i] + shift;
 
1098
    idx    = a->j + a->i[i] ;
 
1099
    v      = a->a + a->i[i] ;
1122
1100
    n      = a->i[i+1] - a->i[i];
1123
1101
    alpha1 = x[8*i];
1124
1102
    alpha2 = x[8*i+1];
1146
1124
  PetscFunctionReturn(0);
1147
1125
}
1148
1126
 
 
1127
/* ------------------------------------------------------------------------------*/
 
1128
#undef __FUNCT__  
 
1129
#define __FUNCT__ "MatMult_SeqMAIJ_9"
 
1130
int MatMult_SeqMAIJ_9(Mat A,Vec xx,Vec yy)
 
1131
{
 
1132
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1133
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1134
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9;
 
1135
  int           ierr,m = b->AIJ->m,*idx,*ii;
 
1136
  int           n,i,jrow,j;
 
1137
 
 
1138
  PetscFunctionBegin;
 
1139
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1140
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
 
1141
  idx  = a->j;
 
1142
  v    = a->a;
 
1143
  ii   = a->i;
 
1144
 
 
1145
  for (i=0; i<m; i++) {
 
1146
    jrow = ii[i];
 
1147
    n    = ii[i+1] - jrow;
 
1148
    sum1  = 0.0;
 
1149
    sum2  = 0.0;
 
1150
    sum3  = 0.0;
 
1151
    sum4  = 0.0;
 
1152
    sum5  = 0.0;
 
1153
    sum6  = 0.0;
 
1154
    sum7  = 0.0;
 
1155
    sum8  = 0.0;
 
1156
    sum9  = 0.0;
 
1157
    for (j=0; j<n; j++) {
 
1158
      sum1 += v[jrow]*x[9*idx[jrow]];
 
1159
      sum2 += v[jrow]*x[9*idx[jrow]+1];
 
1160
      sum3 += v[jrow]*x[9*idx[jrow]+2];
 
1161
      sum4 += v[jrow]*x[9*idx[jrow]+3];
 
1162
      sum5 += v[jrow]*x[9*idx[jrow]+4];
 
1163
      sum6 += v[jrow]*x[9*idx[jrow]+5];
 
1164
      sum7 += v[jrow]*x[9*idx[jrow]+6];
 
1165
      sum8 += v[jrow]*x[9*idx[jrow]+7];
 
1166
      sum9 += v[jrow]*x[9*idx[jrow]+8];
 
1167
      jrow++;
 
1168
     }
 
1169
    y[9*i]   = sum1;
 
1170
    y[9*i+1] = sum2;
 
1171
    y[9*i+2] = sum3;
 
1172
    y[9*i+3] = sum4;
 
1173
    y[9*i+4] = sum5;
 
1174
    y[9*i+5] = sum6;
 
1175
    y[9*i+6] = sum7;
 
1176
    y[9*i+7] = sum8;
 
1177
    y[9*i+8] = sum9;
 
1178
  }
 
1179
 
 
1180
  PetscLogFlops(18*a->nz - 9*m);
 
1181
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1182
  ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
 
1183
  PetscFunctionReturn(0);
 
1184
}
 
1185
 
 
1186
#undef __FUNCT__  
 
1187
#define __FUNCT__ "MatMultTranspose_SeqMAIJ_9"
 
1188
int MatMultTranspose_SeqMAIJ_9(Mat A,Vec xx,Vec yy)
 
1189
{
 
1190
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1191
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1192
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8,alpha9,zero = 0.0;
 
1193
  int           ierr,m = b->AIJ->m,n,i,*idx;
 
1194
 
 
1195
  PetscFunctionBegin; 
 
1196
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
 
1197
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1198
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
 
1199
 
 
1200
  for (i=0; i<m; i++) {
 
1201
    idx    = a->j + a->i[i] ;
 
1202
    v      = a->a + a->i[i] ;
 
1203
    n      = a->i[i+1] - a->i[i];
 
1204
    alpha1 = x[9*i];
 
1205
    alpha2 = x[9*i+1];
 
1206
    alpha3 = x[9*i+2];
 
1207
    alpha4 = x[9*i+3];
 
1208
    alpha5 = x[9*i+4];
 
1209
    alpha6 = x[9*i+5];
 
1210
    alpha7 = x[9*i+6];
 
1211
    alpha8 = x[9*i+7];
 
1212
    alpha9 = x[9*i+8];
 
1213
    while (n-->0) {
 
1214
      y[9*(*idx)]   += alpha1*(*v);
 
1215
      y[9*(*idx)+1] += alpha2*(*v);
 
1216
      y[9*(*idx)+2] += alpha3*(*v);
 
1217
      y[9*(*idx)+3] += alpha4*(*v);
 
1218
      y[9*(*idx)+4] += alpha5*(*v);
 
1219
      y[9*(*idx)+5] += alpha6*(*v);
 
1220
      y[9*(*idx)+6] += alpha7*(*v);
 
1221
      y[9*(*idx)+7] += alpha8*(*v);
 
1222
      y[9*(*idx)+8] += alpha9*(*v);
 
1223
      idx++; v++;
 
1224
    }
 
1225
  }
 
1226
  PetscLogFlops(18*a->nz - 9*b->AIJ->n);
 
1227
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1228
  ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
 
1229
  PetscFunctionReturn(0);
 
1230
}
 
1231
 
 
1232
#undef __FUNCT__  
 
1233
#define __FUNCT__ "MatMultAdd_SeqMAIJ_9"
 
1234
int MatMultAdd_SeqMAIJ_9(Mat A,Vec xx,Vec yy,Vec zz)
 
1235
{
 
1236
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1237
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1238
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9;
 
1239
  int           ierr,m = b->AIJ->m,*idx,*ii;
 
1240
  int           n,i,jrow,j;
 
1241
 
 
1242
  PetscFunctionBegin;
 
1243
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
 
1244
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1245
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
 
1246
  idx  = a->j;
 
1247
  v    = a->a;
 
1248
  ii   = a->i;
 
1249
 
 
1250
  for (i=0; i<m; i++) {
 
1251
    jrow = ii[i];
 
1252
    n    = ii[i+1] - jrow;
 
1253
    sum1  = 0.0;
 
1254
    sum2  = 0.0;
 
1255
    sum3  = 0.0;
 
1256
    sum4  = 0.0;
 
1257
    sum5  = 0.0;
 
1258
    sum6  = 0.0;
 
1259
    sum7  = 0.0;
 
1260
    sum8  = 0.0;
 
1261
    sum9  = 0.0;
 
1262
    for (j=0; j<n; j++) {
 
1263
      sum1 += v[jrow]*x[9*idx[jrow]];
 
1264
      sum2 += v[jrow]*x[9*idx[jrow]+1];
 
1265
      sum3 += v[jrow]*x[9*idx[jrow]+2];
 
1266
      sum4 += v[jrow]*x[9*idx[jrow]+3];
 
1267
      sum5 += v[jrow]*x[9*idx[jrow]+4];
 
1268
      sum6 += v[jrow]*x[9*idx[jrow]+5];
 
1269
      sum7 += v[jrow]*x[9*idx[jrow]+6];
 
1270
      sum8 += v[jrow]*x[9*idx[jrow]+7];
 
1271
      sum9 += v[jrow]*x[9*idx[jrow]+8];
 
1272
      jrow++;
 
1273
     }
 
1274
    y[9*i]   += sum1;
 
1275
    y[9*i+1] += sum2;
 
1276
    y[9*i+2] += sum3;
 
1277
    y[9*i+3] += sum4;
 
1278
    y[9*i+4] += sum5;
 
1279
    y[9*i+5] += sum6;
 
1280
    y[9*i+6] += sum7;
 
1281
    y[9*i+7] += sum8;
 
1282
    y[9*i+8] += sum9;
 
1283
  }
 
1284
 
 
1285
  PetscLogFlops(18*a->nz);
 
1286
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1287
  ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr);
 
1288
  PetscFunctionReturn(0);
 
1289
}
 
1290
 
 
1291
#undef __FUNCT__  
 
1292
#define __FUNCT__ "MatMultTransposeAdd_SeqMAIJ_9"
 
1293
int MatMultTransposeAdd_SeqMAIJ_9(Mat A,Vec xx,Vec yy,Vec zz)
 
1294
{
 
1295
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1296
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1297
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8,alpha9;
 
1298
  int           ierr,m = b->AIJ->m,n,i,*idx;
 
1299
 
 
1300
  PetscFunctionBegin; 
 
1301
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
 
1302
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1303
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
 
1304
  for (i=0; i<m; i++) {
 
1305
    idx    = a->j + a->i[i] ;
 
1306
    v      = a->a + a->i[i] ;
 
1307
    n      = a->i[i+1] - a->i[i];
 
1308
    alpha1 = x[9*i];
 
1309
    alpha2 = x[9*i+1];
 
1310
    alpha3 = x[9*i+2];
 
1311
    alpha4 = x[9*i+3];
 
1312
    alpha5 = x[9*i+4];
 
1313
    alpha6 = x[9*i+5];
 
1314
    alpha7 = x[9*i+6];
 
1315
    alpha8 = x[9*i+7];
 
1316
    alpha9 = x[9*i+8];
 
1317
    while (n-->0) {
 
1318
      y[9*(*idx)]   += alpha1*(*v);
 
1319
      y[9*(*idx)+1] += alpha2*(*v);
 
1320
      y[9*(*idx)+2] += alpha3*(*v);
 
1321
      y[9*(*idx)+3] += alpha4*(*v);
 
1322
      y[9*(*idx)+4] += alpha5*(*v);
 
1323
      y[9*(*idx)+5] += alpha6*(*v);
 
1324
      y[9*(*idx)+6] += alpha7*(*v);
 
1325
      y[9*(*idx)+7] += alpha8*(*v);
 
1326
      y[9*(*idx)+8] += alpha9*(*v);
 
1327
      idx++; v++;
 
1328
    }
 
1329
  }
 
1330
  PetscLogFlops(18*a->nz);
 
1331
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1332
  ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr);
 
1333
  PetscFunctionReturn(0);
 
1334
}
 
1335
 
 
1336
/*--------------------------------------------------------------------------------------------*/
 
1337
#undef __FUNCT__  
 
1338
#define __FUNCT__ "MatMult_SeqMAIJ_16"
 
1339
int MatMult_SeqMAIJ_16(Mat A,Vec xx,Vec yy)
 
1340
{
 
1341
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1342
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1343
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8;
 
1344
  PetscScalar   sum9, sum10, sum11, sum12, sum13, sum14, sum15, sum16;
 
1345
  int           ierr,m = b->AIJ->m,*idx,*ii;
 
1346
  int           n,i,jrow,j;
 
1347
 
 
1348
  PetscFunctionBegin;
 
1349
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1350
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
 
1351
  idx  = a->j;
 
1352
  v    = a->a;
 
1353
  ii   = a->i;
 
1354
 
 
1355
  for (i=0; i<m; i++) {
 
1356
    jrow = ii[i];
 
1357
    n    = ii[i+1] - jrow;
 
1358
    sum1  = 0.0;
 
1359
    sum2  = 0.0;
 
1360
    sum3  = 0.0;
 
1361
    sum4  = 0.0;
 
1362
    sum5  = 0.0;
 
1363
    sum6  = 0.0;
 
1364
    sum7  = 0.0;
 
1365
    sum8  = 0.0;
 
1366
    sum9  = 0.0;
 
1367
    sum10 = 0.0;
 
1368
    sum11 = 0.0;
 
1369
    sum12 = 0.0;
 
1370
    sum13 = 0.0;
 
1371
    sum14 = 0.0;
 
1372
    sum15 = 0.0;
 
1373
    sum16 = 0.0;
 
1374
    for (j=0; j<n; j++) {
 
1375
      sum1  += v[jrow]*x[16*idx[jrow]];
 
1376
      sum2  += v[jrow]*x[16*idx[jrow]+1];
 
1377
      sum3  += v[jrow]*x[16*idx[jrow]+2];
 
1378
      sum4  += v[jrow]*x[16*idx[jrow]+3];
 
1379
      sum5  += v[jrow]*x[16*idx[jrow]+4];
 
1380
      sum6  += v[jrow]*x[16*idx[jrow]+5];
 
1381
      sum7  += v[jrow]*x[16*idx[jrow]+6];
 
1382
      sum8  += v[jrow]*x[16*idx[jrow]+7];
 
1383
      sum9  += v[jrow]*x[16*idx[jrow]+8];
 
1384
      sum10 += v[jrow]*x[16*idx[jrow]+9];
 
1385
      sum11 += v[jrow]*x[16*idx[jrow]+10];
 
1386
      sum12 += v[jrow]*x[16*idx[jrow]+11];
 
1387
      sum13 += v[jrow]*x[16*idx[jrow]+12];
 
1388
      sum14 += v[jrow]*x[16*idx[jrow]+13];
 
1389
      sum15 += v[jrow]*x[16*idx[jrow]+14];
 
1390
      sum16 += v[jrow]*x[16*idx[jrow]+15];
 
1391
      jrow++;
 
1392
     }
 
1393
    y[16*i]    = sum1;
 
1394
    y[16*i+1]  = sum2;
 
1395
    y[16*i+2]  = sum3;
 
1396
    y[16*i+3]  = sum4;
 
1397
    y[16*i+4]  = sum5;
 
1398
    y[16*i+5]  = sum6;
 
1399
    y[16*i+6]  = sum7;
 
1400
    y[16*i+7]  = sum8;
 
1401
    y[16*i+8]  = sum9;
 
1402
    y[16*i+9]  = sum10;
 
1403
    y[16*i+10] = sum11;
 
1404
    y[16*i+11] = sum12;
 
1405
    y[16*i+12] = sum13;
 
1406
    y[16*i+13] = sum14;
 
1407
    y[16*i+14] = sum15;
 
1408
    y[16*i+15] = sum16;
 
1409
  }
 
1410
 
 
1411
  PetscLogFlops(32*a->nz - 16*m);
 
1412
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1413
  ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
 
1414
  PetscFunctionReturn(0);
 
1415
}
 
1416
 
 
1417
#undef __FUNCT__  
 
1418
#define __FUNCT__ "MatMultTranspose_SeqMAIJ_16"
 
1419
int MatMultTranspose_SeqMAIJ_16(Mat A,Vec xx,Vec yy)
 
1420
{
 
1421
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1422
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1423
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8,zero = 0.0;
 
1424
  PetscScalar   alpha9,alpha10,alpha11,alpha12,alpha13,alpha14,alpha15,alpha16;
 
1425
  int           ierr,m = b->AIJ->m,n,i,*idx;
 
1426
 
 
1427
  PetscFunctionBegin; 
 
1428
  ierr = VecSet(&zero,yy);CHKERRQ(ierr);
 
1429
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1430
  ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
 
1431
 
 
1432
  for (i=0; i<m; i++) {
 
1433
    idx    = a->j + a->i[i] ;
 
1434
    v      = a->a + a->i[i] ;
 
1435
    n      = a->i[i+1] - a->i[i];
 
1436
    alpha1  = x[16*i];
 
1437
    alpha2  = x[16*i+1];
 
1438
    alpha3  = x[16*i+2];
 
1439
    alpha4  = x[16*i+3];
 
1440
    alpha5  = x[16*i+4];
 
1441
    alpha6  = x[16*i+5];
 
1442
    alpha7  = x[16*i+6];
 
1443
    alpha8  = x[16*i+7];
 
1444
    alpha9  = x[16*i+8];
 
1445
    alpha10 = x[16*i+9];
 
1446
    alpha11 = x[16*i+10];
 
1447
    alpha12 = x[16*i+11];
 
1448
    alpha13 = x[16*i+12];
 
1449
    alpha14 = x[16*i+13];
 
1450
    alpha15 = x[16*i+14];
 
1451
    alpha16 = x[16*i+15];
 
1452
    while (n-->0) {
 
1453
      y[16*(*idx)]    += alpha1*(*v);
 
1454
      y[16*(*idx)+1]  += alpha2*(*v);
 
1455
      y[16*(*idx)+2]  += alpha3*(*v);
 
1456
      y[16*(*idx)+3]  += alpha4*(*v);
 
1457
      y[16*(*idx)+4]  += alpha5*(*v);
 
1458
      y[16*(*idx)+5]  += alpha6*(*v);
 
1459
      y[16*(*idx)+6]  += alpha7*(*v);
 
1460
      y[16*(*idx)+7]  += alpha8*(*v);
 
1461
      y[16*(*idx)+8]  += alpha9*(*v);
 
1462
      y[16*(*idx)+9]  += alpha10*(*v);
 
1463
      y[16*(*idx)+10] += alpha11*(*v);
 
1464
      y[16*(*idx)+11] += alpha12*(*v);
 
1465
      y[16*(*idx)+12] += alpha13*(*v);
 
1466
      y[16*(*idx)+13] += alpha14*(*v);
 
1467
      y[16*(*idx)+14] += alpha15*(*v);
 
1468
      y[16*(*idx)+15] += alpha16*(*v);
 
1469
      idx++; v++;
 
1470
    }
 
1471
  }
 
1472
  PetscLogFlops(32*a->nz - 16*b->AIJ->n);
 
1473
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1474
  ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
 
1475
  PetscFunctionReturn(0);
 
1476
}
 
1477
 
 
1478
#undef __FUNCT__  
 
1479
#define __FUNCT__ "MatMultAdd_SeqMAIJ_16"
 
1480
int MatMultAdd_SeqMAIJ_16(Mat A,Vec xx,Vec yy,Vec zz)
 
1481
{
 
1482
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1483
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1484
  PetscScalar   *x,*y,*v,sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8;
 
1485
  PetscScalar   sum9, sum10, sum11, sum12, sum13, sum14, sum15, sum16;
 
1486
  int           ierr,m = b->AIJ->m,*idx,*ii;
 
1487
  int           n,i,jrow,j;
 
1488
 
 
1489
  PetscFunctionBegin;
 
1490
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
 
1491
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1492
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
 
1493
  idx  = a->j;
 
1494
  v    = a->a;
 
1495
  ii   = a->i;
 
1496
 
 
1497
  for (i=0; i<m; i++) {
 
1498
    jrow = ii[i];
 
1499
    n    = ii[i+1] - jrow;
 
1500
    sum1  = 0.0;
 
1501
    sum2  = 0.0;
 
1502
    sum3  = 0.0;
 
1503
    sum4  = 0.0;
 
1504
    sum5  = 0.0;
 
1505
    sum6  = 0.0;
 
1506
    sum7  = 0.0;
 
1507
    sum8  = 0.0;
 
1508
    sum9  = 0.0;
 
1509
    sum10 = 0.0;
 
1510
    sum11 = 0.0;
 
1511
    sum12 = 0.0;
 
1512
    sum13 = 0.0;
 
1513
    sum14 = 0.0;
 
1514
    sum15 = 0.0;
 
1515
    sum16 = 0.0;
 
1516
    for (j=0; j<n; j++) {
 
1517
      sum1  += v[jrow]*x[16*idx[jrow]];
 
1518
      sum2  += v[jrow]*x[16*idx[jrow]+1];
 
1519
      sum3  += v[jrow]*x[16*idx[jrow]+2];
 
1520
      sum4  += v[jrow]*x[16*idx[jrow]+3];
 
1521
      sum5  += v[jrow]*x[16*idx[jrow]+4];
 
1522
      sum6  += v[jrow]*x[16*idx[jrow]+5];
 
1523
      sum7  += v[jrow]*x[16*idx[jrow]+6];
 
1524
      sum8  += v[jrow]*x[16*idx[jrow]+7];
 
1525
      sum9  += v[jrow]*x[16*idx[jrow]+8];
 
1526
      sum10 += v[jrow]*x[16*idx[jrow]+9];
 
1527
      sum11 += v[jrow]*x[16*idx[jrow]+10];
 
1528
      sum12 += v[jrow]*x[16*idx[jrow]+11];
 
1529
      sum13 += v[jrow]*x[16*idx[jrow]+12];
 
1530
      sum14 += v[jrow]*x[16*idx[jrow]+13];
 
1531
      sum15 += v[jrow]*x[16*idx[jrow]+14];
 
1532
      sum16 += v[jrow]*x[16*idx[jrow]+15];
 
1533
      jrow++;
 
1534
     }
 
1535
    y[16*i]    += sum1;
 
1536
    y[16*i+1]  += sum2;
 
1537
    y[16*i+2]  += sum3;
 
1538
    y[16*i+3]  += sum4;
 
1539
    y[16*i+4]  += sum5;
 
1540
    y[16*i+5]  += sum6;
 
1541
    y[16*i+6]  += sum7;
 
1542
    y[16*i+7]  += sum8;
 
1543
    y[16*i+8]  += sum9;
 
1544
    y[16*i+9]  += sum10;
 
1545
    y[16*i+10] += sum11;
 
1546
    y[16*i+11] += sum12;
 
1547
    y[16*i+12] += sum13;
 
1548
    y[16*i+13] += sum14;
 
1549
    y[16*i+14] += sum15;
 
1550
    y[16*i+15] += sum16;
 
1551
  }
 
1552
 
 
1553
  PetscLogFlops(32*a->nz);
 
1554
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1555
  ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr);
 
1556
  PetscFunctionReturn(0);
 
1557
}
 
1558
 
 
1559
#undef __FUNCT__  
 
1560
#define __FUNCT__ "MatMultTransposeAdd_SeqMAIJ_16"
 
1561
int MatMultTransposeAdd_SeqMAIJ_16(Mat A,Vec xx,Vec yy,Vec zz)
 
1562
{
 
1563
  Mat_SeqMAIJ   *b = (Mat_SeqMAIJ*)A->data;
 
1564
  Mat_SeqAIJ    *a = (Mat_SeqAIJ*)b->AIJ->data;
 
1565
  PetscScalar   *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,alpha7,alpha8;
 
1566
  PetscScalar   alpha9,alpha10,alpha11,alpha12,alpha13,alpha14,alpha15,alpha16;
 
1567
  int           ierr,m = b->AIJ->m,n,i,*idx;
 
1568
 
 
1569
  PetscFunctionBegin; 
 
1570
  if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);}
 
1571
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
 
1572
  ierr = VecGetArray(zz,&y);CHKERRQ(ierr);
 
1573
  for (i=0; i<m; i++) {
 
1574
    idx    = a->j + a->i[i] ;
 
1575
    v      = a->a + a->i[i] ;
 
1576
    n      = a->i[i+1] - a->i[i];
 
1577
    alpha1 = x[16*i];
 
1578
    alpha2 = x[16*i+1];
 
1579
    alpha3 = x[16*i+2];
 
1580
    alpha4 = x[16*i+3];
 
1581
    alpha5 = x[16*i+4];
 
1582
    alpha6 = x[16*i+5];
 
1583
    alpha7 = x[16*i+6];
 
1584
    alpha8 = x[16*i+7];
 
1585
    alpha9  = x[16*i+8];
 
1586
    alpha10 = x[16*i+9];
 
1587
    alpha11 = x[16*i+10];
 
1588
    alpha12 = x[16*i+11];
 
1589
    alpha13 = x[16*i+12];
 
1590
    alpha14 = x[16*i+13];
 
1591
    alpha15 = x[16*i+14];
 
1592
    alpha16 = x[16*i+15];
 
1593
    while (n-->0) {
 
1594
      y[16*(*idx)]   += alpha1*(*v);
 
1595
      y[16*(*idx)+1] += alpha2*(*v);
 
1596
      y[16*(*idx)+2] += alpha3*(*v);
 
1597
      y[16*(*idx)+3] += alpha4*(*v);
 
1598
      y[16*(*idx)+4] += alpha5*(*v);
 
1599
      y[16*(*idx)+5] += alpha6*(*v);
 
1600
      y[16*(*idx)+6] += alpha7*(*v);
 
1601
      y[16*(*idx)+7] += alpha8*(*v);
 
1602
      y[16*(*idx)+8]  += alpha9*(*v);
 
1603
      y[16*(*idx)+9]  += alpha10*(*v);
 
1604
      y[16*(*idx)+10] += alpha11*(*v);
 
1605
      y[16*(*idx)+11] += alpha12*(*v);
 
1606
      y[16*(*idx)+12] += alpha13*(*v);
 
1607
      y[16*(*idx)+13] += alpha14*(*v);
 
1608
      y[16*(*idx)+14] += alpha15*(*v);
 
1609
      y[16*(*idx)+15] += alpha16*(*v);
 
1610
      idx++; v++;
 
1611
    }
 
1612
  }
 
1613
  PetscLogFlops(32*a->nz);
 
1614
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
 
1615
  ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr);
 
1616
  PetscFunctionReturn(0);
 
1617
}
 
1618
 
1149
1619
/*===================================================================================*/
1150
1620
#undef __FUNCT__  
1151
1621
#define __FUNCT__ "MatMult_MPIMAIJ_dof"
1208
1678
}
1209
1679
 
1210
1680
/* ---------------------------------------------------------------------------------- */
 
1681
/*MC
 
1682
  MatCreateMAIJ - Creates a matrix type providing restriction and interpolation 
 
1683
  operations for multicomponent problems.  It interpolates each component the same
 
1684
  way independently.  The matrix type is based on MATSEQAIJ for sequential matrices,
 
1685
  and MATMPIAIJ for distributed matrices.
 
1686
 
 
1687
  Operations provided:
 
1688
. MatMult
 
1689
. MatMultTranspose
 
1690
. MatMultAdd
 
1691
. MatMultTransposeAdd
 
1692
 
 
1693
  Level: advanced
 
1694
 
 
1695
M*/
1211
1696
#undef __FUNCT__  
1212
1697
#define __FUNCT__ "MatCreateMAIJ" 
1213
1698
int MatCreateMAIJ(Mat A,int dof,Mat *maij)
1262
1747
        B->ops->multadd          = MatMultAdd_SeqMAIJ_8;
1263
1748
        B->ops->multtranspose    = MatMultTranspose_SeqMAIJ_8;
1264
1749
        B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_8;
 
1750
      } else if (dof == 9) {
 
1751
        B->ops->mult             = MatMult_SeqMAIJ_9;
 
1752
        B->ops->multadd          = MatMultAdd_SeqMAIJ_9;
 
1753
        B->ops->multtranspose    = MatMultTranspose_SeqMAIJ_9;
 
1754
        B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_9;
 
1755
      } else if (dof == 16) {
 
1756
        B->ops->mult             = MatMult_SeqMAIJ_16;
 
1757
        B->ops->multadd          = MatMultAdd_SeqMAIJ_16;
 
1758
        B->ops->multtranspose    = MatMultTranspose_SeqMAIJ_16;
 
1759
        B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_16;
1265
1760
      } else {
1266
1761
        SETERRQ1(1,"Cannot handle a dof of %d. Send request for code to petsc-maint@mcs.anl.gov\n",dof);
1267
1762
      }