1 /*$Id: maij.c,v 1.2 2000/05/27 03:52:46 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 17 #include "src/mat/impls/aij/seq/aij.h" 18 19 typedef struct { 20 int dof; /* number of components */ 21 Mat AIJ; /* representation of interpolation for one component */ 22 } Mat_SeqMAIJ; 23 24 #undef __FUNC__ 25 #define __FUNC__ /*<a name="MatDestroy_SeqMAIJ"></a>*/"MatDestroy_SeqMAIJ" 26 int MatDestroy_SeqMAIJ(Mat A) 27 { 28 int ierr; 29 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 30 31 PetscFunctionBegin; 32 if (b->AIJ) { 33 ierr = MatDestroy(b->AIJ);CHKERRQ(ierr); 34 } 35 ierr = PetscFree(b);CHKERRQ(ierr); 36 PetscHeaderDestroy(A); 37 PetscFunctionReturn(0); 38 } 39 40 EXTERN_C_BEGIN 41 #undef __FUNC__ 42 #define __FUNC__ /*<a name="MatCreate_SeqMAIJ"></a>*/"MatCreate_SeqMAIJ" 43 int MatCreate_SeqMAIJ(Mat A) 44 { 45 int ierr; 46 Mat_SeqMAIJ *b; 47 48 PetscFunctionBegin; 49 A->data = (void*)(b = PetscNew(Mat_SeqMAIJ));CHKPTRQ(b); 50 ierr = PetscMemzero(b,sizeof(Mat_SeqMAIJ));CHKERRQ(ierr); 51 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 52 A->factor = 0; 53 A->mapping = 0; 54 b->AIJ = 0; 55 b->dof = 0; 56 PetscFunctionReturn(0); 57 } 58 EXTERN_C_END 59 60 /* ----------------------------------------------------------------------------------*/ 61 EXTERN int MatMult_SeqAIJ(Mat,Vec,Vec); 62 EXTERN int MatMultTranspose_SeqAIJ(Mat,Vec,Vec); 63 EXTERN int MatMultTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 64 EXTERN int matmulttransposeadd_seqaijMatMultAdd_SeqAIJ(Mat,Vec,Vec,Vec); 65 66 #undef __FUNC__ 67 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_1"></a>*/"MatMult_SeqMAIJ_1" 68 int MatMult_SeqMAIJ_1(Mat A,Vec xx,Vec yy) 69 { 70 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 71 int ierr; 72 PetscFunctionBegin; 73 ierr = MatMult_SeqAIJ(b->AIJ,xx,yy); 74 PetscFunctionReturn(0); 75 } 76 #undef __FUNC__ 77 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_1"></a>*/"MatMultTranspose_SeqMAIJ_1" 78 int MatMultTranspose_SeqMAIJ_1(Mat A,Vec xx,Vec yy) 79 { 80 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 81 int ierr; 82 PetscFunctionBegin; 83 ierr = MatMultTranspose_SeqAIJ(b->AIJ,xx,yy); 84 PetscFunctionReturn(0); 85 } 86 #undef __FUNC__ 87 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_1"></a>*/"MatMultAdd_SeqMAIJ_1" 88 int MatMultAdd_SeqMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 89 { 90 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 91 int ierr; 92 PetscFunctionBegin; 93 ierr = MatMultAdd_SeqAIJ(b->AIJ,xx,yy,zz); 94 PetscFunctionReturn(0); 95 } 96 #undef __FUNC__ 97 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_1"></a>*/"MatMultTransposeAdd_SeqMAIJ_1" 98 int MatMultTransposeAdd_SeqMAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 99 { 100 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 101 int ierr; 102 PetscFunctionBegin; 103 ierr = MatMultTransposeAdd_SeqAIJ(b->AIJ,xx,yy,zz); 104 PetscFunctionReturn(0); 105 } 106 107 /* --------------------------------------------------------------------------------------*/ 108 #undef __FUNC__ 109 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_2"></a>*/"MatMult_SeqMAIJ_2" 110 int MatMult_SeqMAIJ_2(Mat A,Vec xx,Vec yy) 111 { 112 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 113 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 114 Scalar *x,*y,*v,sum1, sum2; 115 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 116 int n,i,jrow,j; 117 118 PetscFunctionBegin; 119 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 120 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 121 x = x + shift; /* shift for Fortran start by 1 indexing */ 122 idx = a->j; 123 v = a->a; 124 ii = a->i; 125 126 v += shift; /* shift for Fortran start by 1 indexing */ 127 idx += shift; 128 for (i=0; i<m; i++) { 129 jrow = ii[i]; 130 n = ii[i+1] - jrow; 131 sum1 = 0.0; 132 sum2 = 0.0; 133 for (j=0; j<n; j++) { 134 sum1 += v[jrow]*x[2*idx[jrow]]; 135 sum2 += v[jrow]*x[2*idx[jrow]+1]; 136 jrow++; 137 } 138 y[2*i] = sum1; 139 y[2*i+1] = sum2; 140 } 141 142 PLogFlops(4*a->nz - 2*m); 143 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 144 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNC__ 149 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_2"></a>*/"MatMultTranspose_SeqMAIJ_2" 150 int MatMultTranspose_SeqMAIJ_2(Mat A,Vec xx,Vec yy) 151 { 152 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 153 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 154 Scalar *x,*y,*v,alpha1,alpha2,zero = 0.0; 155 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 156 157 PetscFunctionBegin; 158 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 159 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 160 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 161 y = y + shift; /* shift for Fortran start by 1 indexing */ 162 for (i=0; i<m; i++) { 163 idx = a->j + a->i[i] + shift; 164 v = a->a + a->i[i] + shift; 165 n = a->i[i+1] - a->i[i]; 166 alpha1 = x[2*i]; 167 alpha2 = x[2*i+1]; 168 while (n-->0) {y[2*(*idx)] += alpha1*(*v); y[2*(*idx)+1] += alpha2*(*v); idx++; v++;} 169 } 170 PLogFlops(4*a->nz - 2*a->n); 171 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 172 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNC__ 177 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_2"></a>*/"MatMultAdd_SeqMAIJ_2" 178 int MatMultAdd_SeqMAIJ_2(Mat A,Vec xx,Vec yy,Vec zz) 179 { 180 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 181 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 182 Scalar *x,*y,*v,sum1, sum2; 183 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 184 int n,i,jrow,j; 185 186 PetscFunctionBegin; 187 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 188 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 189 x = x + shift; /* shift for Fortran start by 1 indexing */ 190 idx = a->j; 191 v = a->a; 192 ii = a->i; 193 194 v += shift; /* shift for Fortran start by 1 indexing */ 195 idx += shift; 196 for (i=0; i<m; i++) { 197 jrow = ii[i]; 198 n = ii[i+1] - jrow; 199 sum1 = 0.0; 200 sum2 = 0.0; 201 for (j=0; j<n; j++) { 202 sum1 += v[jrow]*x[2*idx[jrow]]; 203 sum2 += v[jrow]*x[2*idx[jrow]+1]; 204 jrow++; 205 } 206 y[2*i] += sum1; 207 y[2*i+1] += sum2; 208 } 209 210 PLogFlops(4*a->nz - 2*m); 211 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 212 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 PetscFunctionReturn(0); 215 } 216 #undef __FUNC__ 217 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_2"></a>*/"MatMultTransposeAdd_SeqMAIJ_2" 218 int MatMultTransposeAdd_SeqMAIJ_2(Mat A,Vec xx,Vec yy,Vec zz) 219 { 220 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 221 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 222 Scalar *x,*y,*v,alpha1,alpha2; 223 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 224 225 PetscFunctionBegin; 226 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 227 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 228 y = y + shift; /* shift for Fortran start by 1 indexing */ 229 for (i=0; i<m; i++) { 230 idx = a->j + a->i[i] + shift; 231 v = a->a + a->i[i] + shift; 232 n = a->i[i+1] - a->i[i]; 233 alpha1 = x[2*i]; 234 alpha2 = x[2*i+1]; 235 while (n-->0) {y[2*(*idx)] += alpha1*(*v); y[2*(*idx)+1] += alpha2*(*v); idx++; v++;} 236 } 237 PLogFlops(4*a->nz - 2*a->n); 238 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 239 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 PetscFunctionReturn(0); 242 } 243 /* --------------------------------------------------------------------------------------*/ 244 #undef __FUNC__ 245 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_3"></a>*/"MatMult_SeqMAIJ_3" 246 int MatMult_SeqMAIJ_3(Mat A,Vec xx,Vec yy) 247 { 248 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 249 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 250 Scalar *x,*y,*v,sum1, sum2, sum3; 251 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 252 int n,i,jrow,j; 253 254 PetscFunctionBegin; 255 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 256 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 257 x = x + shift; /* shift for Fortran start by 1 indexing */ 258 idx = a->j; 259 v = a->a; 260 ii = a->i; 261 262 v += shift; /* shift for Fortran start by 1 indexing */ 263 idx += shift; 264 for (i=0; i<m; i++) { 265 jrow = ii[i]; 266 n = ii[i+1] - jrow; 267 sum1 = 0.0; 268 sum2 = 0.0; 269 sum3 = 0.0; 270 for (j=0; j<n; j++) { 271 sum1 += v[jrow]*x[3*idx[jrow]]; 272 sum2 += v[jrow]*x[3*idx[jrow]+1]; 273 sum3 += v[jrow]*x[3*idx[jrow]+2]; 274 jrow++; 275 } 276 y[3*i] = sum1; 277 y[3*i+1] = sum2; 278 y[3*i+2] = sum3; 279 } 280 281 PLogFlops(6*a->nz - 3*m); 282 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 283 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 284 PetscFunctionReturn(0); 285 } 286 287 #undef __FUNC__ 288 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_3"></a>*/"MatMultTranspose_SeqMAIJ_3" 289 int MatMultTranspose_SeqMAIJ_3(Mat A,Vec xx,Vec yy) 290 { 291 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 292 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 293 Scalar *x,*y,*v,alpha1,alpha2,alpha3,zero = 0.0; 294 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 295 296 PetscFunctionBegin; 297 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 298 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 299 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 300 y = y + shift; /* shift for Fortran start by 1 indexing */ 301 for (i=0; i<m; i++) { 302 idx = a->j + a->i[i] + shift; 303 v = a->a + a->i[i] + shift; 304 n = a->i[i+1] - a->i[i]; 305 alpha1 = x[3*i]; 306 alpha2 = x[3*i+1]; 307 alpha3 = x[3*i+2]; 308 while (n-->0) { 309 y[3*(*idx)] += alpha1*(*v); 310 y[3*(*idx)+1] += alpha2*(*v); 311 y[3*(*idx)+2] += alpha3*(*v); 312 idx++; v++; 313 } 314 } 315 PLogFlops(6*a->nz - 3*a->n); 316 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 317 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 } 320 321 #undef __FUNC__ 322 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_3"></a>*/"MatMultAdd_SeqMAIJ_3" 323 int MatMultAdd_SeqMAIJ_3(Mat A,Vec xx,Vec yy,Vec zz) 324 { 325 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 326 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 327 Scalar *x,*y,*v,sum1, sum2, sum3; 328 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 329 int n,i,jrow,j; 330 331 PetscFunctionBegin; 332 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 333 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 334 x = x + shift; /* shift for Fortran start by 1 indexing */ 335 idx = a->j; 336 v = a->a; 337 ii = a->i; 338 339 v += shift; /* shift for Fortran start by 1 indexing */ 340 idx += shift; 341 for (i=0; i<m; i++) { 342 jrow = ii[i]; 343 n = ii[i+1] - jrow; 344 sum1 = 0.0; 345 sum2 = 0.0; 346 sum3 = 0.0; 347 for (j=0; j<n; j++) { 348 sum1 += v[jrow]*x[3*idx[jrow]]; 349 sum2 += v[jrow]*x[3*idx[jrow]+1]; 350 sum3 += v[jrow]*x[3*idx[jrow]+2]; 351 jrow++; 352 } 353 y[3*i] += sum1; 354 y[3*i+1] += sum2; 355 y[3*i+2] += sum3; 356 } 357 358 PLogFlops(6*a->nz); 359 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 360 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 #undef __FUNC__ 364 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_3"></a>*/"MatMultTransposeAdd_SeqMAIJ_3" 365 int MatMultTransposeAdd_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,alpha1,alpha2,alpha3; 370 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 371 372 PetscFunctionBegin; 373 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 374 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 375 y = y + shift; /* shift for Fortran start by 1 indexing */ 376 for (i=0; i<m; i++) { 377 idx = a->j + a->i[i] + shift; 378 v = a->a + a->i[i] + shift; 379 n = a->i[i+1] - a->i[i]; 380 alpha1 = x[3*i]; 381 alpha2 = x[3*i+1]; 382 alpha3 = x[3*i+2]; 383 while (n-->0) { 384 y[3*(*idx)] += alpha1*(*v); 385 y[3*(*idx)+1] += alpha2*(*v); 386 y[3*(*idx)+2] += alpha3*(*v); 387 idx++; v++; 388 } 389 } 390 PLogFlops(6*a->nz); 391 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 392 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 /* ------------------------------------------------------------------------------*/ 397 #undef __FUNC__ 398 #define __FUNC__ /*<a name="MatMult_SeqMAIJ_4"></a>*/"MatMult_SeqMAIJ_4" 399 int MatMult_SeqMAIJ_4(Mat A,Vec xx,Vec yy) 400 { 401 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 402 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 403 Scalar *x,*y,*v,sum1, sum2, sum3, sum4; 404 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 405 int n,i,jrow,j; 406 407 PetscFunctionBegin; 408 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 409 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 410 x = x + shift; /* shift for Fortran start by 1 indexing */ 411 idx = a->j; 412 v = a->a; 413 ii = a->i; 414 415 v += shift; /* shift for Fortran start by 1 indexing */ 416 idx += shift; 417 for (i=0; i<m; i++) { 418 jrow = ii[i]; 419 n = ii[i+1] - jrow; 420 sum1 = 0.0; 421 sum2 = 0.0; 422 sum3 = 0.0; 423 sum4 = 0.0; 424 for (j=0; j<n; j++) { 425 sum1 += v[jrow]*x[4*idx[jrow]]; 426 sum2 += v[jrow]*x[4*idx[jrow]+1]; 427 sum3 += v[jrow]*x[4*idx[jrow]+2]; 428 sum4 += v[jrow]*x[4*idx[jrow]+3]; 429 jrow++; 430 } 431 y[4*i] = sum1; 432 y[4*i+1] = sum2; 433 y[4*i+2] = sum3; 434 y[4*i+3] = sum4; 435 } 436 437 PLogFlops(8*a->nz - 4*m); 438 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 439 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 443 #undef __FUNC__ 444 #define __FUNC__ /*<a name="MatMultTranspose_SeqMAIJ_4"></a>*/"MatMultTranspose_SeqMAIJ_4" 445 int MatMultTranspose_SeqMAIJ_4(Mat A,Vec xx,Vec yy) 446 { 447 Mat_SeqMAIJ *b = (Mat_SeqMAIJ*)A->data; 448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)b->AIJ->data; 449 Scalar *x,*y,*v,alpha1,alpha2,alpha3,alpha4,zero = 0.0; 450 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 451 452 PetscFunctionBegin; 453 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 454 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 455 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 456 y = y + shift; /* shift for Fortran start by 1 indexing */ 457 for (i=0; i<m; i++) { 458 idx = a->j + a->i[i] + shift; 459 v = a->a + a->i[i] + shift; 460 n = a->i[i+1] - a->i[i]; 461 alpha1 = x[4*i]; 462 alpha2 = x[4*i+1]; 463 alpha3 = x[4*i+2]; 464 alpha4 = x[4*i+3]; 465 while (n-->0) { 466 y[4*(*idx)] += alpha1*(*v); 467 y[4*(*idx)+1] += alpha2*(*v); 468 y[4*(*idx)+2] += alpha3*(*v); 469 y[4*(*idx)+3] += alpha4*(*v); 470 idx++; v++; 471 } 472 } 473 PLogFlops(8*a->nz - 4*a->n); 474 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 475 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 476 PetscFunctionReturn(0); 477 } 478 479 #undef __FUNC__ 480 #define __FUNC__ /*<a name="MatMultAdd_SeqMAIJ_4"></a>*/"MatMultAdd_SeqMAIJ_4" 481 int MatMultAdd_SeqMAIJ_4(Mat A,Vec xx,Vec yy,Vec zz) 482 { 483 return 0; 484 } 485 #undef __FUNC__ 486 #define __FUNC__ /*<a name="MatMultTransposeAdd_SeqMAIJ_4"></a>*/"MatMultTransposeAdd_SeqMAIJ_4" 487 int MatMultTransposeAdd_SeqMAIJ_4(Mat A,Vec xx,Vec yy,Vec zz) 488 { 489 return 0; 490 } 491 492 /* ---------------------------------------------------------------------------------- */ 493 #undef __FUNC__ 494 #define __FUNC__ /*<a name="MatCreateMAIJ"></a>*/"MatCreateMAIJ" 495 int MatCreateMAIJ(Mat A,int dof,Mat *maij) 496 { 497 int ierr; 498 Mat_SeqMAIJ *b; 499 Mat B; 500 501 PetscFunctionBegin; 502 ierr = MATCreate(A->comm,dof*A->m,dof*A->n,dof*A->M,dof*A->N,&B);CHKERRQ(ierr); 503 ierr = MatSetType(B,MATSEQMAIJ);CHKERRQ(ierr); 504 505 B->assembled = PETSC_TRUE; 506 B->ops->destroy = MatDestroy_SeqMAIJ; 507 b = (Mat_SeqMAIJ*)B->data; 508 509 b->AIJ = A; 510 b->dof = dof; 511 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 512 if (dof == 1) { 513 B->ops->mult = MatMult_SeqMAIJ_1; 514 B->ops->multadd = MatMultAdd_SeqMAIJ_1; 515 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_1; 516 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_1; 517 } else if (dof == 2) { 518 B->ops->mult = MatMult_SeqMAIJ_2; 519 B->ops->multadd = MatMultAdd_SeqMAIJ_2; 520 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_2; 521 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_2; 522 } else if (dof == 3) { 523 B->ops->mult = MatMult_SeqMAIJ_3; 524 B->ops->multadd = MatMultAdd_SeqMAIJ_3; 525 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_3; 526 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_3; 527 } else if (dof == 4) { 528 B->ops->mult = MatMult_SeqMAIJ_4; 529 B->ops->multadd = MatMultAdd_SeqMAIJ_4; 530 B->ops->multtranspose = MatMultTranspose_SeqMAIJ_4; 531 B->ops->multtransposeadd = MatMultTransposeAdd_SeqMAIJ_4; 532 } else { 533 SETERRQ1(1,1,"Cannot handle a dof of %d\n",dof); 534 } 535 *maij = B; 536 PetscFunctionReturn(0); 537 } 538 539 540 541 542 543 544 545 546 547 548 549 550