1 /*$Id: maij.c,v 1.5 2000/06/05 19:38:31 bsmith Exp bsmith $*/ 2 /* 3 Defines the basic matrix operations for the MAIJ matrix storage format. 4 This format is used for restriction and interpolation operations for 5 multicomponent problems. It interpolates each component the same way 6 independently. 7 8 We provide: 9 MatMult() 10 MatMultTranspose() 11 MatMultTransposeAdd() 12 MatMultAdd() 13 and 14 MatCreateMAIJ(Mat,dof,Mat*) 15 16 This single directory handles both the sequential and parallel codes 17 */ 18 19 #include "src/mat/impls/aij/mpi/mpiaij.h" 20 21 typedef struct { 22 int dof; /* number of components */ 23 Mat AIJ; /* representation of interpolation for one component */ 24 } Mat_SeqMAIJ; 25 26 typedef struct { 27 int dof; /* number of components */ 28 Mat AIJ,OAIJ; /* representation of interpolation for one component */ 29 Mat A; 30 VecScatter ctx; /* update ghost points for parallel case */ 31 Vec w; /* work space for ghost values for parallel case */ 32 } Mat_MPIMAIJ; 33 34 #undef __FUNC__ 35 #define __FUNC__ /*<a name="MatDestroy_SeqMAIJ"></a>*/"MatDestroy_SeqMAIJ" 36 int MatDestroy_SeqMAIJ(Mat A) 37 { 38 int ierr; 39 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 40 41 PetscFunctionBegin; 42 if (b->AIJ) { 43 ierr = MatDestroy(b->AIJ);CHKERRQ(ierr); 44 } 45 ierr = PetscFree(b);CHKERRQ(ierr); 46 PetscHeaderDestroy(A); 47 PetscFunctionReturn(0); 48 } 49 50 #undef __FUNC__ 51 #define __FUNC__ /*<a name="MatDestroy_MPIMAIJ"></a>*/"MatDestroy_MPIMAIJ" 52 int MatDestroy_MPIMAIJ(Mat A) 53 { 54 int ierr; 55 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 56 57 PetscFunctionBegin; 58 if (b->A) { 59 ierr = MatDestroy(b->A);CHKERRQ(ierr); 60 } 61 if (b->AIJ) { 62 ierr = MatDestroy(b->AIJ);CHKERRQ(ierr); 63 } 64 if (b->OAIJ) { 65 ierr = MatDestroy(b->OAIJ);CHKERRQ(ierr); 66 } 67 if (b->ctx) { 68 ierr = VecScatterDestroy(b->ctx);CHKERRQ(ierr); 69 } 70 if (b->w) { 71 ierr = VecDestroy(b->w);CHKERRQ(ierr); 72 } 73 ierr = PetscFree(b);CHKERRQ(ierr); 74 PetscHeaderDestroy(A); 75 PetscFunctionReturn(0); 76 } 77 78 EXTERN_C_BEGIN 79 #undef __FUNC__ 80 #define __FUNC__ /*<a name="MatCreate_MAIJ"></a>*/"MatCreate_MAIJ" 81 int MatCreate_MAIJ(Mat A) 82 { 83 int ierr; 84 Mat_MPIMAIJ *b; 85 86 PetscFunctionBegin; 87 A->data = (void*)(b = PetscNew(Mat_MPIMAIJ));CHKPTRQ(b); 88 ierr = PetscMemzero(b,sizeof(Mat_MPIMAIJ));CHKERRQ(ierr); 89 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 90 A->factor = 0; 91 A->mapping = 0; 92 93 b->AIJ = 0; 94 b->dof = 0; 95 b->OAIJ = 0; 96 b->ctx = 0; 97 b->w = 0; 98 PetscFunctionReturn(0); 99 } 100 EXTERN_C_END 101 102 /* ----------------------------------------------------------------------------------*/ 103 EXTERN int MatMult_SeqAIJ(Mat,Vec,Vec); 104 EXTERN int MatMultAdd_SeqAIJ(Mat,Vec,Vec,Vec); 105 EXTERN int MatMultTranspose_SeqAIJ(Mat,Vec,Vec); 106 EXTERN int MatMultTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 107 108 #undef __FUNC__ 109 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_1"></a>*/"MatMult_SeqMAIJ_1" 110 int MatMult_SeqMAIJ_1(Mat A,Vec xx,Vec yy) 111 { 112 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 113 int ierr; 114 PetscFunctionBegin; 115 ierr = MatMult_SeqAIJ(b->AIJ,xx,yy); 116 PetscFunctionReturn(0); 117 } 118 #undef __FUNC__ 119 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_1"></a>*/"MatMultTranspose_SeqMAIJ_1" 120 int MatMultTranspose_SeqMAIJ_1(Mat A,Vec xx,Vec yy) 121 { 122 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 123 int ierr; 124 PetscFunctionBegin; 125 ierr = MatMultTranspose_SeqAIJ(b->AIJ,xx,yy); 126 PetscFunctionReturn(0); 127 } 128 #undef __FUNC__ 129 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_1"></a>*/"MatMultAdd_SeqMAIJ_1" 130 int MatMultAdd_SeqMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 131 { 132 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 133 int ierr; 134 PetscFunctionBegin; 135 ierr = MatMultAdd_SeqAIJ(b->AIJ,xx,yy,zz); 136 PetscFunctionReturn(0); 137 } 138 #undef __FUNC__ 139 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_1"></a>*/"MatMultTransposeAdd_SeqMAIJ_1" 140 int MatMultTransposeAdd_SeqMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 141 { 142 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 143 int ierr; 144 PetscFunctionBegin; 145 ierr = MatMultTransposeAdd_SeqAIJ(b->AIJ,xx,yy,zz); 146 PetscFunctionReturn(0); 147 } 148 149 /* --------------------------------------------------------------------------------------*/ 150 #undef __FUNC__ 151 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_2"></a>*/"MatMult_SeqMAIJ_2" 152 int MatMult_SeqMAIJ_2(Mat A,Vec xx,Vec yy) 153 { 154 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 155 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 156 Scalar *x,*y,*v,sum1, sum2; 157 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 158 int n,i,jrow,j; 159 160 PetscFunctionBegin; 161 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 162 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 163 x = x + shift; /* shift for Fortran start by 1 indexing */ 164 idx = a->j; 165 v = a->a; 166 ii = a->i; 167 168 v += shift; /* shift for Fortran start by 1 indexing */ 169 idx += shift; 170 for (i=0; i<m; i++) { 171 jrow = ii[i]; 172 n = ii[i+1] - jrow; 173 sum1 = 0.0; 174 sum2 = 0.0; 175 for (j=0; j<n; j++) { 176 sum1 += v[jrow]*x[2*idx[jrow]]; 177 sum2 += v[jrow]*x[2*idx[jrow]+1]; 178 jrow++; 179 } 180 y[2*i] = sum1; 181 y[2*i+1] = sum2; 182 } 183 184 PLogFlops(4*a->nz - 2*m); 185 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 186 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNC__ 191 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_2"></a>*/"MatMultTranspose_SeqMAIJ_2" 192 int MatMultTranspose_SeqMAIJ_2(Mat A,Vec xx,Vec yy) 193 { 194 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 195 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 196 Scalar *x,*y,*v,alpha1,alpha2,zero = 0.0; 197 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 198 199 PetscFunctionBegin; 200 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 201 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 202 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 203 y = y + shift; /* shift for Fortran start by 1 indexing */ 204 for (i=0; i<m; i++) { 205 idx = a->j + a->i[i] + shift; 206 v = a->a + a->i[i] + shift; 207 n = a->i[i+1] - a->i[i]; 208 alpha1 = x[2*i]; 209 alpha2 = x[2*i+1]; 210 while (n-->0) {y[2*(*idx)] += alpha1*(*v); y[2*(*idx)+1] += alpha2*(*v); idx++; v++;} 211 } 212 PLogFlops(4*a->nz - 2*a->n); 213 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 214 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNC__ 219 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_2"></a>*/"MatMultAdd_SeqMAIJ_2" 220 int MatMultAdd_SeqMAIJ_2(Mat A,Vec xx,Vec yy,Vec zz) 221 { 222 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 224 Scalar *x,*y,*v,sum1, sum2; 225 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 226 int n,i,jrow,j; 227 228 PetscFunctionBegin; 229 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 230 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 231 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 232 x = x + shift; /* shift for Fortran start by 1 indexing */ 233 idx = a->j; 234 v = a->a; 235 ii = a->i; 236 237 v += shift; /* shift for Fortran start by 1 indexing */ 238 idx += shift; 239 for (i=0; i<m; i++) { 240 jrow = ii[i]; 241 n = ii[i+1] - jrow; 242 sum1 = 0.0; 243 sum2 = 0.0; 244 for (j=0; j<n; j++) { 245 sum1 += v[jrow]*x[2*idx[jrow]]; 246 sum2 += v[jrow]*x[2*idx[jrow]+1]; 247 jrow++; 248 } 249 y[2*i] += sum1; 250 y[2*i+1] += sum2; 251 } 252 253 PLogFlops(4*a->nz - 2*m); 254 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 255 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 #undef __FUNC__ 259 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_2"></a>*/"MatMultTransposeAdd_SeqMAIJ_2" 260 int MatMultTransposeAdd_SeqMAIJ_2(Mat A,Vec xx,Vec yy,Vec zz) 261 { 262 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 263 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 264 Scalar *x,*y,*v,alpha1,alpha2; 265 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 266 267 PetscFunctionBegin; 268 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 269 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 270 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 271 y = y + shift; /* shift for Fortran start by 1 indexing */ 272 for (i=0; i<m; i++) { 273 idx = a->j + a->i[i] + shift; 274 v = a->a + a->i[i] + shift; 275 n = a->i[i+1] - a->i[i]; 276 alpha1 = x[2*i]; 277 alpha2 = x[2*i+1]; 278 while (n-->0) {y[2*(*idx)] += alpha1*(*v); y[2*(*idx)+1] += alpha2*(*v); idx++; v++;} 279 } 280 PLogFlops(4*a->nz - 2*a->n); 281 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 282 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 /* --------------------------------------------------------------------------------------*/ 286 #undef __FUNC__ 287 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_3"></a>*/"MatMult_SeqMAIJ_3" 288 int MatMult_SeqMAIJ_3(Mat A,Vec xx,Vec yy) 289 { 290 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 291 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 292 Scalar *x,*y,*v,sum1, sum2, sum3; 293 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 294 int n,i,jrow,j; 295 296 PetscFunctionBegin; 297 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 298 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 299 x = x + shift; /* shift for Fortran start by 1 indexing */ 300 idx = a->j; 301 v = a->a; 302 ii = a->i; 303 304 v += shift; /* shift for Fortran start by 1 indexing */ 305 idx += shift; 306 for (i=0; i<m; i++) { 307 jrow = ii[i]; 308 n = ii[i+1] - jrow; 309 sum1 = 0.0; 310 sum2 = 0.0; 311 sum3 = 0.0; 312 for (j=0; j<n; j++) { 313 sum1 += v[jrow]*x[3*idx[jrow]]; 314 sum2 += v[jrow]*x[3*idx[jrow]+1]; 315 sum3 += v[jrow]*x[3*idx[jrow]+2]; 316 jrow++; 317 } 318 y[3*i] = sum1; 319 y[3*i+1] = sum2; 320 y[3*i+2] = sum3; 321 } 322 323 PLogFlops(6*a->nz - 3*m); 324 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 325 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328 329 #undef __FUNC__ 330 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_3"></a>*/"MatMultTranspose_SeqMAIJ_3" 331 int MatMultTranspose_SeqMAIJ_3(Mat A,Vec xx,Vec yy) 332 { 333 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 334 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 335 Scalar *x,*y,*v,alpha1,alpha2,alpha3,zero = 0.0; 336 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 337 338 PetscFunctionBegin; 339 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 340 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 341 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 342 y = y + shift; /* shift for Fortran start by 1 indexing */ 343 for (i=0; i<m; i++) { 344 idx = a->j + a->i[i] + shift; 345 v = a->a + a->i[i] + shift; 346 n = a->i[i+1] - a->i[i]; 347 alpha1 = x[3*i]; 348 alpha2 = x[3*i+1]; 349 alpha3 = x[3*i+2]; 350 while (n-->0) { 351 y[3*(*idx)] += alpha1*(*v); 352 y[3*(*idx)+1] += alpha2*(*v); 353 y[3*(*idx)+2] += alpha3*(*v); 354 idx++; v++; 355 } 356 } 357 PLogFlops(6*a->nz - 3*a->n); 358 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 359 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 #undef __FUNC__ 364 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_3"></a>*/"MatMultAdd_SeqMAIJ_3" 365 int MatMultAdd_SeqMAIJ_3(Mat A,Vec xx,Vec yy,Vec zz) 366 { 367 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 368 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 369 Scalar *x,*y,*v,sum1, sum2, sum3; 370 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 371 int n,i,jrow,j; 372 373 PetscFunctionBegin; 374 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 375 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 376 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 377 x = x + shift; /* shift for Fortran start by 1 indexing */ 378 idx = a->j; 379 v = a->a; 380 ii = a->i; 381 382 v += shift; /* shift for Fortran start by 1 indexing */ 383 idx += shift; 384 for (i=0; i<m; i++) { 385 jrow = ii[i]; 386 n = ii[i+1] - jrow; 387 sum1 = 0.0; 388 sum2 = 0.0; 389 sum3 = 0.0; 390 for (j=0; j<n; j++) { 391 sum1 += v[jrow]*x[3*idx[jrow]]; 392 sum2 += v[jrow]*x[3*idx[jrow]+1]; 393 sum3 += v[jrow]*x[3*idx[jrow]+2]; 394 jrow++; 395 } 396 y[3*i] += sum1; 397 y[3*i+1] += sum2; 398 y[3*i+2] += sum3; 399 } 400 401 PLogFlops(6*a->nz); 402 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 403 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 404 PetscFunctionReturn(0); 405 } 406 #undef __FUNC__ 407 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_3"></a>*/"MatMultTransposeAdd_SeqMAIJ_3" 408 int MatMultTransposeAdd_SeqMAIJ_3(Mat A,Vec xx,Vec yy,Vec zz) 409 { 410 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 411 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 412 Scalar *x,*y,*v,alpha1,alpha2,alpha3; 413 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 414 415 PetscFunctionBegin; 416 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 417 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 418 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 419 y = y + shift; /* shift for Fortran start by 1 indexing */ 420 for (i=0; i<m; i++) { 421 idx = a->j + a->i[i] + shift; 422 v = a->a + a->i[i] + shift; 423 n = a->i[i+1] - a->i[i]; 424 alpha1 = x[3*i]; 425 alpha2 = x[3*i+1]; 426 alpha3 = x[3*i+2]; 427 while (n-->0) { 428 y[3*(*idx)] += alpha1*(*v); 429 y[3*(*idx)+1] += alpha2*(*v); 430 y[3*(*idx)+2] += alpha3*(*v); 431 idx++; v++; 432 } 433 } 434 PLogFlops(6*a->nz); 435 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 436 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 437 PetscFunctionReturn(0); 438 } 439 440 /* ------------------------------------------------------------------------------*/ 441 #undef __FUNC__ 442 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_4"></a>*/"MatMult_SeqMAIJ_4" 443 int MatMult_SeqMAIJ_4(Mat A,Vec xx,Vec yy) 444 { 445 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 446 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 447 Scalar *x,*y,*v,sum1, sum2, sum3, sum4; 448 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 449 int n,i,jrow,j; 450 451 PetscFunctionBegin; 452 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 453 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 454 x = x + shift; /* shift for Fortran start by 1 indexing */ 455 idx = a->j; 456 v = a->a; 457 ii = a->i; 458 459 v += shift; /* shift for Fortran start by 1 indexing */ 460 idx += shift; 461 for (i=0; i<m; i++) { 462 jrow = ii[i]; 463 n = ii[i+1] - jrow; 464 sum1 = 0.0; 465 sum2 = 0.0; 466 sum3 = 0.0; 467 sum4 = 0.0; 468 for (j=0; j<n; j++) { 469 sum1 += v[jrow]*x[4*idx[jrow]]; 470 sum2 += v[jrow]*x[4*idx[jrow]+1]; 471 sum3 += v[jrow]*x[4*idx[jrow]+2]; 472 sum4 += v[jrow]*x[4*idx[jrow]+3]; 473 jrow++; 474 } 475 y[4*i] = sum1; 476 y[4*i+1] = sum2; 477 y[4*i+2] = sum3; 478 y[4*i+3] = sum4; 479 } 480 481 PLogFlops(8*a->nz - 4*m); 482 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 483 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 484 PetscFunctionReturn(0); 485 } 486 487 #undef __FUNC__ 488 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_4"></a>*/"MatMultTranspose_SeqMAIJ_4" 489 int MatMultTranspose_SeqMAIJ_4(Mat A,Vec xx,Vec yy) 490 { 491 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 492 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 493 Scalar *x,*y,*v,alpha1,alpha2,alpha3,alpha4,zero = 0.0; 494 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 495 496 PetscFunctionBegin; 497 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 498 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 499 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 500 y = y + shift; /* shift for Fortran start by 1 indexing */ 501 for (i=0; i<m; i++) { 502 idx = a->j + a->i[i] + shift; 503 v = a->a + a->i[i] + shift; 504 n = a->i[i+1] - a->i[i]; 505 alpha1 = x[4*i]; 506 alpha2 = x[4*i+1]; 507 alpha3 = x[4*i+2]; 508 alpha4 = x[4*i+3]; 509 while (n-->0) { 510 y[4*(*idx)] += alpha1*(*v); 511 y[4*(*idx)+1] += alpha2*(*v); 512 y[4*(*idx)+2] += alpha3*(*v); 513 y[4*(*idx)+3] += alpha4*(*v); 514 idx++; v++; 515 } 516 } 517 PLogFlops(8*a->nz - 4*a->n); 518 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 519 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 523 #undef __FUNC__ 524 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_4"></a>*/"MatMultAdd_SeqMAIJ_4" 525 int MatMultAdd_SeqMAIJ_4(Mat A,Vec xx,Vec yy,Vec zz) 526 { 527 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 528 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 529 Scalar *x,*y,*v,sum1, sum2, sum3, sum4; 530 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 531 int n,i,jrow,j; 532 533 PetscFunctionBegin; 534 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 535 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 536 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 537 x = x + shift; /* shift for Fortran start by 1 indexing */ 538 idx = a->j; 539 v = a->a; 540 ii = a->i; 541 542 v += shift; /* shift for Fortran start by 1 indexing */ 543 idx += shift; 544 for (i=0; i<m; i++) { 545 jrow = ii[i]; 546 n = ii[i+1] - jrow; 547 sum1 = 0.0; 548 sum2 = 0.0; 549 sum3 = 0.0; 550 sum4 = 0.0; 551 for (j=0; j<n; j++) { 552 sum1 += v[jrow]*x[4*idx[jrow]]; 553 sum2 += v[jrow]*x[4*idx[jrow]+1]; 554 sum3 += v[jrow]*x[4*idx[jrow]+2]; 555 sum4 += v[jrow]*x[4*idx[jrow]+3]; 556 jrow++; 557 } 558 y[4*i] += sum1; 559 y[4*i+1] += sum2; 560 y[4*i+2] += sum3; 561 y[4*i+3] += sum4; 562 } 563 564 PLogFlops(8*a->nz - 4*m); 565 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 566 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 567 PetscFunctionReturn(0); 568 } 569 #undef __FUNC__ 570 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_4"></a>*/"MatMultTransposeAdd_SeqMAIJ_4" 571 int MatMultTransposeAdd_SeqMAIJ_4(Mat A,Vec xx,Vec yy,Vec zz) 572 { 573 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 574 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 575 Scalar *x,*y,*v,alpha1,alpha2,alpha3,alpha4; 576 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 577 578 PetscFunctionBegin; 579 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 580 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 581 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 582 y = y + shift; /* shift for Fortran start by 1 indexing */ 583 for (i=0; i<m; i++) { 584 idx = a->j + a->i[i] + shift; 585 v = a->a + a->i[i] + shift; 586 n = a->i[i+1] - a->i[i]; 587 alpha1 = x[4*i]; 588 alpha2 = x[4*i+1]; 589 alpha3 = x[4*i+2]; 590 alpha4 = x[4*i+3]; 591 while (n-->0) { 592 y[4*(*idx)] += alpha1*(*v); 593 y[4*(*idx)+1] += alpha2*(*v); 594 y[4*(*idx)+2] += alpha3*(*v); 595 y[4*(*idx)+3] += alpha4*(*v); 596 idx++; v++; 597 } 598 } 599 PLogFlops(8*a->nz - 4*a->n); 600 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 601 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 602 PetscFunctionReturn(0); 603 604 } 605 606 /* ------------------------------------------------------------------------------*/ 607 #undef __FUNC__ 608 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_5"></a>*/"MatMult_SeqMAIJ_5" 609 int MatMult_SeqMAIJ_5(Mat A,Vec xx,Vec yy) 610 { 611 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 612 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 613 Scalar *x,*y,*v,sum1, sum2, sum3, sum4, sum5; 614 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 615 int n,i,jrow,j; 616 617 PetscFunctionBegin; 618 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 619 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 620 x = x + shift; /* shift for Fortran start by 1 indexing */ 621 idx = a->j; 622 v = a->a; 623 ii = a->i; 624 625 v += shift; /* shift for Fortran start by 1 indexing */ 626 idx += shift; 627 for (i=0; i<m; i++) { 628 jrow = ii[i]; 629 n = ii[i+1] - jrow; 630 sum1 = 0.0; 631 sum2 = 0.0; 632 sum3 = 0.0; 633 sum4 = 0.0; 634 sum5 = 0.0; 635 for (j=0; j<n; j++) { 636 sum1 += v[jrow]*x[5*idx[jrow]]; 637 sum2 += v[jrow]*x[5*idx[jrow]+1]; 638 sum3 += v[jrow]*x[5*idx[jrow]+2]; 639 sum4 += v[jrow]*x[5*idx[jrow]+3]; 640 sum5 += v[jrow]*x[5*idx[jrow]+4]; 641 jrow++; 642 } 643 y[5*i] = sum1; 644 y[5*i+1] = sum2; 645 y[5*i+2] = sum3; 646 y[5*i+3] = sum4; 647 y[5*i+4] = sum5; 648 } 649 650 PLogFlops(10*a->nz - 5*m); 651 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 652 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 #undef __FUNC__ 657 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_5"></a>*/"MatMultTranspose_SeqMAIJ_5" 658 int MatMultTranspose_SeqMAIJ_5(Mat A,Vec xx,Vec yy) 659 { 660 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 661 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 662 Scalar *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5,zero = 0.0; 663 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 664 665 PetscFunctionBegin; 666 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 667 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 668 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 669 y = y + shift; /* shift for Fortran start by 1 indexing */ 670 for (i=0; i<m; i++) { 671 idx = a->j + a->i[i] + shift; 672 v = a->a + a->i[i] + shift; 673 n = a->i[i+1] - a->i[i]; 674 alpha1 = x[5*i]; 675 alpha2 = x[5*i+1]; 676 alpha3 = x[5*i+2]; 677 alpha4 = x[5*i+3]; 678 alpha5 = x[5*i+4]; 679 while (n-->0) { 680 y[5*(*idx)] += alpha1*(*v); 681 y[5*(*idx)+1] += alpha2*(*v); 682 y[5*(*idx)+2] += alpha3*(*v); 683 y[5*(*idx)+3] += alpha4*(*v); 684 y[5*(*idx)+4] += alpha5*(*v); 685 idx++; v++; 686 } 687 } 688 PLogFlops(10*a->nz - 5*a->n); 689 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 690 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 691 PetscFunctionReturn(0); 692 } 693 694 #undef __FUNC__ 695 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_5"></a>*/"MatMultAdd_SeqMAIJ_5" 696 int MatMultAdd_SeqMAIJ_5(Mat A,Vec xx,Vec yy,Vec zz) 697 { 698 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 699 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 700 Scalar *x,*y,*v,sum1, sum2, sum3, sum4, sum5; 701 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 702 int n,i,jrow,j; 703 704 PetscFunctionBegin; 705 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 706 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 707 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 708 x = x + shift; /* shift for Fortran start by 1 indexing */ 709 idx = a->j; 710 v = a->a; 711 ii = a->i; 712 713 v += shift; /* shift for Fortran start by 1 indexing */ 714 idx += shift; 715 for (i=0; i<m; i++) { 716 jrow = ii[i]; 717 n = ii[i+1] - jrow; 718 sum1 = 0.0; 719 sum2 = 0.0; 720 sum3 = 0.0; 721 sum4 = 0.0; 722 sum5 = 0.0; 723 for (j=0; j<n; j++) { 724 sum1 += v[jrow]*x[5*idx[jrow]]; 725 sum2 += v[jrow]*x[5*idx[jrow]+1]; 726 sum3 += v[jrow]*x[5*idx[jrow]+2]; 727 sum4 += v[jrow]*x[5*idx[jrow]+3]; 728 sum5 += v[jrow]*x[5*idx[jrow]+4]; 729 jrow++; 730 } 731 y[5*i] += sum1; 732 y[5*i+1] += sum2; 733 y[5*i+2] += sum3; 734 y[5*i+3] += sum4; 735 y[5*i+4] += sum5; 736 } 737 738 PLogFlops(10*a->nz); 739 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 740 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNC__ 745 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_5"></a>*/"MatMultTransposeAdd_SeqMAIJ_5" 746 int MatMultTransposeAdd_SeqMAIJ_5(Mat A,Vec xx,Vec yy,Vec zz) 747 { 748 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 749 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 750 Scalar *x,*y,*v,alpha1,alpha2,alpha3,alpha4,alpha5; 751 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 752 753 PetscFunctionBegin; 754 if (yy != zz) {ierr = VecCopy(yy,zz);CHKERRQ(ierr);} 755 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 756 ierr = VecGetArray(zz,&y);CHKERRQ(ierr); 757 y = y + shift; /* shift for Fortran start by 1 indexing */ 758 for (i=0; i<m; i++) { 759 idx = a->j + a->i[i] + shift; 760 v = a->a + a->i[i] + shift; 761 n = a->i[i+1] - a->i[i]; 762 alpha1 = x[5*i]; 763 alpha2 = x[5*i+1]; 764 alpha3 = x[5*i+2]; 765 alpha4 = x[5*i+3]; 766 alpha5 = x[5*i+4]; 767 while (n-->0) { 768 y[5*(*idx)] += alpha1*(*v); 769 y[5*(*idx)+1] += alpha2*(*v); 770 y[5*(*idx)+2] += alpha3*(*v); 771 y[5*(*idx)+3] += alpha4*(*v); 772 y[5*(*idx)+4] += alpha5*(*v); 773 idx++; v++; 774 } 775 } 776 PLogFlops(10*a->nz); 777 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 778 ierr = VecRestoreArray(zz,&y);CHKERRQ(ierr); 779 PetscFunctionReturn(0); 780 } 781 782 /*===================================================================================*/ 783 EXTERN int MatMult_MPIAIJ(Mat,Vec,Vec); 784 EXTERN int MatMultAdd_MPIAIJ(Mat,Vec,Vec,Vec); 785 EXTERN int MatMultTranspose_MPIAIJ(Mat,Vec,Vec); 786 EXTERN int MatMultTransposeAdd_MPIAIJ(Mat,Vec,Vec,Vec); 787 788 #undef __FUNC__ 789 #define __FUNC__ /*<a name="MatMult_MPIMAIJ_1"></a>*/"MatMult_MPIMAIJ_1" 790 int MatMult_MPIMAIJ_1(Mat A,Vec xx,Vec yy) 791 { 792 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 793 int ierr; 794 PetscFunctionBegin; 795 ierr = MatMult_MPIAIJ(b->AIJ,xx,yy); 796 PetscFunctionReturn(0); 797 } 798 #undef __FUNC__ 799 #define __FUNC__ /*<a name="MatMultTranspose_MPIMAIJ_1"></a>*/"MatMultTranspose_MPIMAIJ_1" 800 int MatMultTranspose_MPIMAIJ_1(Mat A,Vec xx,Vec yy) 801 { 802 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 803 int ierr; 804 PetscFunctionBegin; 805 ierr = MatMultTranspose_MPIAIJ(b->AIJ,xx,yy); 806 PetscFunctionReturn(0); 807 } 808 #undef __FUNC__ 809 #define __FUNC__ /*<a name="MatMultAdd_MPIMAIJ_1"></a>*/"MatMultAdd_MPIMAIJ_1" 810 int MatMultAdd_MPIMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 811 { 812 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 813 int ierr; 814 PetscFunctionBegin; 815 ierr = MatMultAdd_MPIAIJ(b->AIJ,xx,yy,zz); 816 PetscFunctionReturn(0); 817 } 818 #undef __FUNC__ 819 #define __FUNC__ /*<a name="MatMultTransposeAdd_MPIMAIJ_1"></a>*/"MatMultTransposeAdd_MPIMAIJ_1" 820 int MatMultTransposeAdd_MPIMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 821 { 822 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 823 int ierr; 824 PetscFunctionBegin; 825 ierr = MatMultTransposeAdd_MPIAIJ(b->AIJ,xx,yy,zz); 826 PetscFunctionReturn(0); 827 } 828 829 /* ---------------------------------------------------------------------------------- */ 830 #undef __FUNC__ 831 #define __FUNC__ /*<a name="MatMult_MPIMAIJ_dof"></a>*/"MatMult_MPIMAIJ_dof" 832 int MatMult_MPIMAIJ_dof(Mat A,Vec xx,Vec yy) 833 { 834 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 835 int ierr; 836 PetscFunctionBegin; 837 838 /* start the scatter */ 839 ierr = VecScatterBegin(xx,b->w,INSERT_VALUES,SCATTER_FORWARD,b->ctx);CHKERRQ(ierr); 840 ierr = (*b->AIJ->ops->mult)(b->AIJ,xx,yy);CHKERRQ(ierr); 841 ierr = VecScatterEnd(xx,b->w,INSERT_VALUES,SCATTER_FORWARD,b->ctx);CHKERRQ(ierr); 842 ierr = (*b->OAIJ->ops->multadd)(b->OAIJ,b->w,yy,yy);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 846 #undef __FUNC__ 847 #define __FUNC__ /*<a name="MatMultTranspose_MPIMAIJ_dof"></a>*/"MatMultTranspose_MPIMAIJ_dof" 848 int MatMultTranspose_MPIMAIJ_dof(Mat A,Vec xx,Vec yy) 849 { 850 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 851 int ierr; 852 PetscFunctionBegin; 853 ierr = (*b->OAIJ->ops->multtranspose)(b->OAIJ,xx,b->w);CHKERRQ(ierr); 854 ierr = VecScatterBegin(b->w,yy,ADD_VALUES,SCATTER_REVERSE,b->ctx);CHKERRQ(ierr); 855 ierr = (*b->AIJ->ops->multtranspose)(b->AIJ,xx,yy);CHKERRQ(ierr); 856 ierr = VecScatterEnd(b->w,yy,ADD_VALUES,SCATTER_REVERSE,b->ctx);CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 #undef __FUNC__ 861 #define __FUNC__ /*<a name="MatMultAdd_MPIMAIJ_dof"></a>*/"MatMultAdd_MPIMAIJ_dof" 862 int MatMultAdd_MPIMAIJ_dof(Mat A,Vec xx,Vec yy) 863 { 864 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 865 int ierr; 866 PetscFunctionBegin; 867 868 /* start the scatter */ 869 ierr = VecScatterBegin(xx,b->w,INSERT_VALUES,SCATTER_FORWARD,b->ctx);CHKERRQ(ierr); 870 ierr = (*b->AIJ->ops->multadd)(b->AIJ,xx,yy,yy);CHKERRQ(ierr); 871 ierr = VecScatterEnd(xx,b->w,INSERT_VALUES,SCATTER_FORWARD,b->ctx);CHKERRQ(ierr); 872 ierr = (*b->OAIJ->ops->multadd)(b->OAIJ,b->w,yy,yy);CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 #undef __FUNC__ 877 #define __FUNC__ /*<a name="MatMultTransposeAdd_MPIMAIJ_dof"></a>*/"MatMultTransposeAdd_MPIMAIJ_dof" 878 int MatMultTransposeAdd_MPIMAIJ_dof(Mat A,Vec xx,Vec yy) 879 { 880 Mat_MPIMAIJ *b = (Mat_MPIMAIJ*)A->data; 881 int ierr; 882 PetscFunctionBegin; 883 ierr = (*b->OAIJ->ops->multtranspose)(b->OAIJ,xx,b->w);CHKERRQ(ierr); 884 ierr = VecScatterBegin(b->w,yy,ADD_VALUES,SCATTER_REVERSE,b->ctx);CHKERRQ(ierr); 885 ierr = (*b->AIJ->ops->multtransposeadd)(b->AIJ,xx,yy,yy);CHKERRQ(ierr); 886 ierr = VecScatterEnd(b->w,yy,ADD_VALUES,SCATTER_REVERSE,b->ctx);CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 891 /* ---------------------------------------------------------------------------------- */ 892 #undef __FUNC__ 893 #define __FUNC__ /*<a name="MatCreateMAIJ"></a>*/"MatCreateMAIJ" 894 int MatCreateMAIJ(Mat A,int dof,Mat *maij) 895 { 896 int ierr,size,n; 897 Mat_MPIMAIJ *b; 898 Mat B; 899 900 PetscFunctionBegin; 901 ierr = MATCreate(A->comm,dof*A->m,dof*A->n,dof*A->M,dof*A->N,&B);CHKERRQ(ierr); 902 ierr = MatSetType(B,MATMAIJ);CHKERRQ(ierr); 903 904 B->assembled = PETSC_TRUE; 905 b = (Mat_MPIMAIJ*)B->data; 906 907 b->dof = dof; 908 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 909 ierr = MPI_Comm_size(A->comm,&size);CHKERRQ(ierr); 910 if (size == 1) { 911 B->ops->destroy = MatDestroy_SeqMAIJ; 912 b->AIJ = A; 913 if (dof == 1) { 914 B->ops->mult = MatMult_SeqMAIJ_1; 915 B->ops->multadd = MatMultAdd_SeqMAIJ_1; 916 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_1; 917 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_1; 918 } else if (dof == 2) { 919 B->ops->mult = MatMult_SeqMAIJ_2; 920 B->ops->multadd = MatMultAdd_SeqMAIJ_2; 921 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_2; 922 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_2; 923 } else if (dof == 3) { 924 B->ops->mult = MatMult_SeqMAIJ_3; 925 B->ops->multadd = MatMultAdd_SeqMAIJ_3; 926 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_3; 927 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_3; 928 } else if (dof == 4) { 929 B->ops->mult = MatMult_SeqMAIJ_4; 930 B->ops->multadd = MatMultAdd_SeqMAIJ_4; 931 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_4; 932 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_4; 933 } else if (dof == 5) { 934 B->ops->mult = MatMult_SeqMAIJ_5; 935 B->ops->multadd = MatMultAdd_SeqMAIJ_5; 936 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_5; 937 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_5; 938 } else { 939 SETERRQ1(1,1,"Cannot handle a dof of %d\n",dof); 940 } 941 } else { 942 b->A = A; 943 B->ops->destroy = MatDestroy_MPIMAIJ; 944 if (dof == 1) { 945 b->AIJ = A; 946 B->ops->mult = MatMult_MPIMAIJ_1; 947 B->ops->multadd = MatMultAdd_MPIMAIJ_1; 948 B->ops->multtranspose = MatMultTranspose_MPIMAIJ_1; 949 B->ops->multtransposeadd = MatMultTransposeAdd_MPIMAIJ_1; 950 } else { 951 Mat_MPIAIJ *mpiaij = (Mat_MPIAIJ *)A->data; 952 IS from,to; 953 Vec gvec; 954 int *garray,i; 955 956 ierr = MatCreateMAIJ(mpiaij->A,dof,&b->AIJ);CHKERRQ(ierr); 957 ierr = MatCreateMAIJ(mpiaij->B,dof,&b->OAIJ);CHKERRQ(ierr); 958 959 ierr = VecGetSize(mpiaij->lvec,&n);CHKERRQ(ierr); 960 ierr = VecCreateSeq(PETSC_COMM_SELF,n*dof,&b->w);CHKERRQ(ierr); 961 962 /* create two temporary Index sets for build scatter gather */ 963 garray = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(garray); 964 for (i=0; i<n; i++) garray[i] = dof*mpiaij->garray[i]; 965 ierr = ISCreateBlock(A->comm,dof,n,garray,&from);CHKERRQ(ierr); 966 ierr = PetscFree(garray);CHKERRQ(ierr); 967 ierr = ISCreateStride(PETSC_COMM_SELF,n*dof,0,1,&to);CHKERRQ(ierr); 968 969 /* create temporary global vector to generate scatter context */ 970 ierr = VecCreateMPI(A->comm,dof*A->n,dof*A->N,&gvec);CHKERRQ(ierr); 971 972 /* generate the scatter context */ 973 ierr = VecScatterCreate(gvec,from,b->w,to,&b->ctx);CHKERRQ(ierr); 974 975 ierr = ISDestroy(from);CHKERRQ(ierr); 976 ierr = ISDestroy(to);CHKERRQ(ierr); 977 ierr = VecDestroy(gvec);CHKERRQ(ierr); 978 979 B->ops->mult = MatMult_MPIMAIJ_dof; 980 B->ops->multtranspose = MatMultTranspose_MPIMAIJ_dof; 981 B->ops->multadd = MatMultAdd_MPIMAIJ_dof; 982 B->ops->multtransposeadd = MatMultTransposeAdd_MPIMAIJ_dof; 983 } 984 } 985 *maij = B; 986 PetscFunctionReturn(0); 987 } 988 989 990 991 992 993 994 995 996 997 998 999 1000