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