1 2 #include "src/mat/impls/baij/seq/baij.h" 3 #include "src/inline/spops.h" 4 #include "src/inline/ilu.h" 5 #include "petscbt.h" 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatIncreaseOverlap_SeqBAIJ" 9 PetscErrorCode MatIncreaseOverlap_SeqBAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 10 { 11 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 12 PetscErrorCode ierr; 13 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val,ival; 14 PetscInt start,end,*ai,*aj,bs,*nidx2; 15 PetscBT table; 16 17 PetscFunctionBegin; 18 m = a->mbs; 19 ai = a->i; 20 aj = a->j; 21 bs = A->bs; 22 23 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified"); 24 25 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 26 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 27 ierr = PetscMalloc((A->m+1)*sizeof(PetscInt),&nidx2);CHKERRQ(ierr); 28 29 for (i=0; i<is_max; i++) { 30 /* Initialise the two local arrays */ 31 isz = 0; 32 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 33 34 /* Extract the indices, assume there can be duplicate entries */ 35 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 36 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 37 38 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 39 for (j=0; j<n ; ++j){ 40 ival = idx[j]/bs; /* convert the indices into block indices */ 41 if (ival>=m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than mat-dim"); 42 if(!PetscBTLookupSet(table,ival)) { nidx[isz++] = ival;} 43 } 44 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 45 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 46 47 k = 0; 48 for (j=0; j<ov; j++){ /* for each overlap*/ 49 n = isz; 50 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 51 row = nidx[k]; 52 start = ai[row]; 53 end = ai[row+1]; 54 for (l = start; l<end ; l++){ 55 val = aj[l]; 56 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 57 } 58 } 59 } 60 /* expand the Index Set */ 61 for (j=0; j<isz; j++) { 62 for (k=0; k<bs; k++) 63 nidx2[j*bs+k] = nidx[j]*bs+k; 64 } 65 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz*bs,nidx2,is+i);CHKERRQ(ierr); 66 } 67 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 68 ierr = PetscFree(nidx);CHKERRQ(ierr); 69 ierr = PetscFree(nidx2);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "MatGetSubMatrix_SeqBAIJ_Private" 75 PetscErrorCode MatGetSubMatrix_SeqBAIJ_Private(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 76 { 77 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*c; 78 PetscErrorCode ierr; 79 PetscInt *smap,i,k,kstart,kend,oldcols = a->nbs,*lens; 80 PetscInt row,mat_i,*mat_j,tcol,*mat_ilen; 81 PetscInt *irow,*icol,nrows,ncols,*ssmap,bs=A->bs,bs2=a->bs2; 82 PetscInt *aj = a->j,*ai = a->i; 83 MatScalar *mat_a; 84 Mat C; 85 PetscTruth flag; 86 87 PetscFunctionBegin; 88 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 89 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IS is not sorted"); 90 91 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 92 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 93 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 94 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 95 96 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 97 ssmap = smap; 98 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 99 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 100 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 101 /* determine lens of each row */ 102 for (i=0; i<nrows; i++) { 103 kstart = ai[irow[i]]; 104 kend = kstart + a->ilen[irow[i]]; 105 lens[i] = 0; 106 for (k=kstart; k<kend; k++) { 107 if (ssmap[aj[k]]) { 108 lens[i]++; 109 } 110 } 111 } 112 /* Create and fill new matrix */ 113 if (scall == MAT_REUSE_MATRIX) { 114 c = (Mat_SeqBAIJ *)((*B)->data); 115 116 if (c->mbs!=nrows || c->nbs!=ncols || (*B)->bs!=bs) SETERRQ(PETSC_ERR_ARG_SIZ,"Submatrix wrong size"); 117 ierr = PetscMemcmp(c->ilen,lens,c->mbs *sizeof(PetscInt),&flag);CHKERRQ(ierr); 118 if (flag == PETSC_FALSE) { 119 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 120 } 121 ierr = PetscMemzero(c->ilen,c->mbs*sizeof(PetscInt));CHKERRQ(ierr); 122 C = *B; 123 } else { 124 ierr = MatCreate(A->comm,nrows*bs,ncols*bs,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 125 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 126 ierr = MatSeqBAIJSetPreallocation(C,bs,0,lens);CHKERRQ(ierr); 127 } 128 c = (Mat_SeqBAIJ *)(C->data); 129 for (i=0; i<nrows; i++) { 130 row = irow[i]; 131 kstart = ai[row]; 132 kend = kstart + a->ilen[row]; 133 mat_i = c->i[i]; 134 mat_j = c->j + mat_i; 135 mat_a = c->a + mat_i*bs2; 136 mat_ilen = c->ilen + i; 137 for (k=kstart; k<kend; k++) { 138 if ((tcol=ssmap[a->j[k]])) { 139 *mat_j++ = tcol - 1; 140 ierr = PetscMemcpy(mat_a,a->a+k*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 141 mat_a += bs2; 142 (*mat_ilen)++; 143 } 144 } 145 } 146 147 /* Free work space */ 148 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 149 ierr = PetscFree(smap);CHKERRQ(ierr); 150 ierr = PetscFree(lens);CHKERRQ(ierr); 151 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 152 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 153 154 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 155 *B = C; 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "MatGetSubMatrix_SeqBAIJ" 161 PetscErrorCode MatGetSubMatrix_SeqBAIJ(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 162 { 163 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 164 IS is1,is2; 165 PetscErrorCode ierr; 166 PetscInt *vary,*iary,*irow,*icol,nrows,ncols,i,bs=A->bs,count; 167 168 PetscFunctionBegin; 169 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 170 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 171 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 172 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 173 174 /* Verify if the indices corespond to each element in a block 175 and form the IS with compressed IS */ 176 ierr = PetscMalloc(2*(a->mbs+1)*sizeof(PetscInt),&vary);CHKERRQ(ierr); 177 iary = vary + a->mbs; 178 ierr = PetscMemzero(vary,(a->mbs)*sizeof(PetscInt));CHKERRQ(ierr); 179 for (i=0; i<nrows; i++) vary[irow[i]/bs]++; 180 count = 0; 181 for (i=0; i<a->mbs; i++) { 182 if (vary[i]!=0 && vary[i]!=bs) SETERRQ(PETSC_ERR_ARG_SIZ,"Index set does not match blocks"); 183 if (vary[i]==bs) iary[count++] = i; 184 } 185 ierr = ISCreateGeneral(PETSC_COMM_SELF,count,iary,&is1);CHKERRQ(ierr); 186 187 ierr = PetscMemzero(vary,(a->mbs)*sizeof(PetscInt));CHKERRQ(ierr); 188 for (i=0; i<ncols; i++) vary[icol[i]/bs]++; 189 count = 0; 190 for (i=0; i<a->mbs; i++) { 191 if (vary[i]!=0 && vary[i]!=bs) SETERRQ(PETSC_ERR_PLIB,"Internal error in PETSc"); 192 if (vary[i]==bs) iary[count++] = i; 193 } 194 ierr = ISCreateGeneral(PETSC_COMM_SELF,count,iary,&is2);CHKERRQ(ierr); 195 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 196 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 197 ierr = PetscFree(vary);CHKERRQ(ierr); 198 199 ierr = MatGetSubMatrix_SeqBAIJ_Private(A,is1,is2,cs,scall,B);CHKERRQ(ierr); 200 ISDestroy(is1); 201 ISDestroy(is2); 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNCT__ 206 #define __FUNCT__ "MatGetSubMatrices_SeqBAIJ" 207 PetscErrorCode MatGetSubMatrices_SeqBAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 208 { 209 PetscErrorCode ierr; 210 PetscInt i; 211 212 PetscFunctionBegin; 213 if (scall == MAT_INITIAL_MATRIX) { 214 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 215 } 216 217 for (i=0; i<n; i++) { 218 ierr = MatGetSubMatrix_SeqBAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 219 } 220 PetscFunctionReturn(0); 221 } 222 223 224 /* -------------------------------------------------------*/ 225 /* Should check that shapes of vectors and matrices match */ 226 /* -------------------------------------------------------*/ 227 #include "petscblaslapack.h" 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatMult_SeqBAIJ_1" 231 PetscErrorCode MatMult_SeqBAIJ_1(Mat A,Vec xx,Vec zz) 232 { 233 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 234 PetscScalar *x,*z,sum; 235 MatScalar *v; 236 PetscErrorCode ierr; 237 PetscInt mbs,i,*idx,*ii,n,*rindex=PETSC_NULL; 238 PetscTruth usecprow=a->compressedrow.use; 239 240 PetscFunctionBegin; 241 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 242 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 243 244 idx = a->j; 245 v = a->a; 246 if (usecprow){ 247 mbs = a->compressedrow.nrows; 248 ii = a->compressedrow.i; 249 rindex = a->compressedrow.rindex; 250 } else { 251 mbs = a->mbs; 252 ii = a->i; 253 } 254 255 for (i=0; i<mbs; i++) { 256 n = ii[1] - ii[0]; ii++; 257 sum = 0.0; 258 while (n--) sum += *v++ * x[*idx++]; 259 if (usecprow){ 260 z[rindex[i]] = sum; 261 } else { 262 z[i] = sum; 263 } 264 } 265 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 266 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 267 PetscLogFlops(2*a->nz - mbs); 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "MatMult_SeqBAIJ_2" 273 PetscErrorCode MatMult_SeqBAIJ_2(Mat A,Vec xx,Vec zz) 274 { 275 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 276 PetscScalar *x,*z,*xb,sum1,sum2,*zarray; 277 PetscScalar x1,x2; 278 MatScalar *v; 279 PetscErrorCode ierr; 280 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 281 PetscTruth usecprow=a->compressedrow.use; 282 283 PetscFunctionBegin; 284 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 285 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 286 287 idx = a->j; 288 v = a->a; 289 if (usecprow){ 290 mbs = a->compressedrow.nrows; 291 ii = a->compressedrow.i; 292 rindex = a->compressedrow.rindex; 293 } else { 294 mbs = a->mbs; 295 ii = a->i; 296 z = zarray; 297 } 298 299 for (i=0; i<mbs; i++) { 300 n = ii[1] - ii[0]; ii++; 301 sum1 = 0.0; sum2 = 0.0; 302 for (j=0; j<n; j++) { 303 xb = x + 2*(*idx++); x1 = xb[0]; x2 = xb[1]; 304 sum1 += v[0]*x1 + v[2]*x2; 305 sum2 += v[1]*x1 + v[3]*x2; 306 v += 4; 307 } 308 if (usecprow) z = zarray + 2*rindex[i]; 309 z[0] = sum1; z[1] = sum2; 310 if (!usecprow) z += 2; 311 } 312 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 313 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 314 PetscLogFlops(8*a->nz - 2*mbs); 315 PetscFunctionReturn(0); 316 } 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatMult_SeqBAIJ_3" 320 PetscErrorCode MatMult_SeqBAIJ_3(Mat A,Vec xx,Vec zz) 321 { 322 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 323 PetscScalar *x,*z,*xb,sum1,sum2,sum3,x1,x2,x3,*zarray; 324 MatScalar *v; 325 PetscErrorCode ierr; 326 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 327 PetscTruth usecprow=a->compressedrow.use; 328 329 330 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 331 #pragma disjoint(*v,*z,*xb) 332 #endif 333 334 PetscFunctionBegin; 335 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 336 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 337 338 idx = a->j; 339 v = a->a; 340 if (usecprow){ 341 mbs = a->compressedrow.nrows; 342 ii = a->compressedrow.i; 343 rindex = a->compressedrow.rindex; 344 } else { 345 mbs = a->mbs; 346 ii = a->i; 347 z = zarray; 348 } 349 350 for (i=0; i<mbs; i++) { 351 n = ii[1] - ii[0]; ii++; 352 sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; 353 for (j=0; j<n; j++) { 354 xb = x + 3*(*idx++); x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 355 sum1 += v[0]*x1 + v[3]*x2 + v[6]*x3; 356 sum2 += v[1]*x1 + v[4]*x2 + v[7]*x3; 357 sum3 += v[2]*x1 + v[5]*x2 + v[8]*x3; 358 v += 9; 359 } 360 if (usecprow) z = zarray + 3*rindex[i]; 361 z[0] = sum1; z[1] = sum2; z[2] = sum3; 362 if (!usecprow) z += 3; 363 } 364 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 365 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 366 PetscLogFlops(18*a->nz - 3*mbs); 367 PetscFunctionReturn(0); 368 } 369 370 #undef __FUNCT__ 371 #define __FUNCT__ "MatMult_SeqBAIJ_4" 372 PetscErrorCode MatMult_SeqBAIJ_4(Mat A,Vec xx,Vec zz) 373 { 374 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 375 PetscScalar *x,*z,*xb,sum1,sum2,sum3,sum4,x1,x2,x3,x4,*zarray; 376 MatScalar *v; 377 PetscErrorCode ierr; 378 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 379 PetscTruth usecprow=a->compressedrow.use; 380 381 PetscFunctionBegin; 382 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 383 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 384 385 idx = a->j; 386 v = a->a; 387 if (usecprow){ 388 mbs = a->compressedrow.nrows; 389 ii = a->compressedrow.i; 390 rindex = a->compressedrow.rindex; 391 } else { 392 mbs = a->mbs; 393 ii = a->i; 394 z = zarray; 395 } 396 397 for (i=0; i<mbs; i++) { 398 n = ii[1] - ii[0]; ii++; 399 sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; sum4 = 0.0; 400 for (j=0; j<n; j++) { 401 xb = x + 4*(*idx++); 402 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; 403 sum1 += v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 404 sum2 += v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 405 sum3 += v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 406 sum4 += v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 407 v += 16; 408 } 409 if (usecprow) z = zarray + 4*rindex[i]; 410 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; 411 if (!usecprow) z += 4; 412 } 413 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 414 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 415 PetscLogFlops(32*a->nz - 4*mbs); 416 PetscFunctionReturn(0); 417 } 418 419 #undef __FUNCT__ 420 #define __FUNCT__ "MatMult_SeqBAIJ_5" 421 PetscErrorCode MatMult_SeqBAIJ_5(Mat A,Vec xx,Vec zz) 422 { 423 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 424 PetscScalar sum1,sum2,sum3,sum4,sum5,x1,x2,x3,x4,x5,*xb,*z,*x,*zarray; 425 MatScalar *v; 426 PetscErrorCode ierr; 427 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 428 PetscTruth usecprow=a->compressedrow.use; 429 430 PetscFunctionBegin; 431 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 432 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 433 434 idx = a->j; 435 v = a->a; 436 if (usecprow){ 437 mbs = a->compressedrow.nrows; 438 ii = a->compressedrow.i; 439 rindex = a->compressedrow.rindex; 440 } else { 441 mbs = a->mbs; 442 ii = a->i; 443 z = zarray; 444 } 445 446 for (i=0; i<mbs; i++) { 447 n = ii[1] - ii[0]; ii++; 448 sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; sum4 = 0.0; sum5 = 0.0; 449 for (j=0; j<n; j++) { 450 xb = x + 5*(*idx++); 451 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; 452 sum1 += v[0]*x1 + v[5]*x2 + v[10]*x3 + v[15]*x4 + v[20]*x5; 453 sum2 += v[1]*x1 + v[6]*x2 + v[11]*x3 + v[16]*x4 + v[21]*x5; 454 sum3 += v[2]*x1 + v[7]*x2 + v[12]*x3 + v[17]*x4 + v[22]*x5; 455 sum4 += v[3]*x1 + v[8]*x2 + v[13]*x3 + v[18]*x4 + v[23]*x5; 456 sum5 += v[4]*x1 + v[9]*x2 + v[14]*x3 + v[19]*x4 + v[24]*x5; 457 v += 25; 458 } 459 if (usecprow) z = zarray + 5*rindex[i]; 460 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; 461 if (!usecprow) z += 5; 462 } 463 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 464 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 465 PetscLogFlops(50*a->nz - 5*mbs); 466 PetscFunctionReturn(0); 467 } 468 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "MatMult_SeqBAIJ_6" 472 PetscErrorCode MatMult_SeqBAIJ_6(Mat A,Vec xx,Vec zz) 473 { 474 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 475 PetscScalar *x,*z,*xb,sum1,sum2,sum3,sum4,sum5,sum6; 476 PetscScalar x1,x2,x3,x4,x5,x6,*zarray; 477 MatScalar *v; 478 PetscErrorCode ierr; 479 PetscInt mbs=a->mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 480 PetscTruth usecprow=a->compressedrow.use; 481 482 PetscFunctionBegin; 483 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 484 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 485 486 idx = a->j; 487 v = a->a; 488 if (usecprow){ 489 mbs = a->compressedrow.nrows; 490 ii = a->compressedrow.i; 491 rindex = a->compressedrow.rindex; 492 } else { 493 mbs = a->mbs; 494 ii = a->i; 495 z = zarray; 496 } 497 498 for (i=0; i<mbs; i++) { 499 n = ii[1] - ii[0]; ii++; 500 sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; sum4 = 0.0; sum5 = 0.0; sum6 = 0.0; 501 for (j=0; j<n; j++) { 502 xb = x + 6*(*idx++); 503 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; 504 sum1 += v[0]*x1 + v[6]*x2 + v[12]*x3 + v[18]*x4 + v[24]*x5 + v[30]*x6; 505 sum2 += v[1]*x1 + v[7]*x2 + v[13]*x3 + v[19]*x4 + v[25]*x5 + v[31]*x6; 506 sum3 += v[2]*x1 + v[8]*x2 + v[14]*x3 + v[20]*x4 + v[26]*x5 + v[32]*x6; 507 sum4 += v[3]*x1 + v[9]*x2 + v[15]*x3 + v[21]*x4 + v[27]*x5 + v[33]*x6; 508 sum5 += v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4 + v[28]*x5 + v[34]*x6; 509 sum6 += v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4 + v[29]*x5 + v[35]*x6; 510 v += 36; 511 } 512 if (usecprow) z = zarray + 6*rindex[i]; 513 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; z[5] = sum6; 514 if (!usecprow) z += 6; 515 } 516 517 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 518 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 519 PetscLogFlops(72*a->nz - 6*mbs); 520 PetscFunctionReturn(0); 521 } 522 #undef __FUNCT__ 523 #define __FUNCT__ "MatMult_SeqBAIJ_7" 524 PetscErrorCode MatMult_SeqBAIJ_7(Mat A,Vec xx,Vec zz) 525 { 526 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 527 PetscScalar *x,*z,*xb,sum1,sum2,sum3,sum4,sum5,sum6,sum7; 528 PetscScalar x1,x2,x3,x4,x5,x6,x7,*zarray; 529 MatScalar *v; 530 PetscErrorCode ierr; 531 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 532 PetscTruth usecprow=a->compressedrow.use; 533 534 PetscFunctionBegin; 535 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 536 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 537 538 idx = a->j; 539 v = a->a; 540 if (usecprow){ 541 mbs = a->compressedrow.nrows; 542 ii = a->compressedrow.i; 543 rindex = a->compressedrow.rindex; 544 } else { 545 mbs = a->mbs; 546 ii = a->i; 547 z = zarray; 548 } 549 550 for (i=0; i<mbs; i++) { 551 n = ii[1] - ii[0]; ii++; 552 sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; sum4 = 0.0; sum5 = 0.0; sum6 = 0.0; sum7 = 0.0; 553 for (j=0; j<n; j++) { 554 xb = x + 7*(*idx++); 555 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; x7 = xb[6]; 556 sum1 += v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 557 sum2 += v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 558 sum3 += v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 559 sum4 += v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 560 sum5 += v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 561 sum6 += v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 562 sum7 += v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 563 v += 49; 564 } 565 if (usecprow) z = zarray + 7*rindex[i]; 566 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; z[5] = sum6; z[6] = sum7; 567 if (!usecprow) z += 7; 568 } 569 570 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 571 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 572 PetscLogFlops(98*a->nz - 7*mbs); 573 PetscFunctionReturn(0); 574 } 575 576 /* 577 This will not work with MatScalar == float because it calls the BLAS 578 */ 579 #undef __FUNCT__ 580 #define __FUNCT__ "MatMult_SeqBAIJ_N" 581 PetscErrorCode MatMult_SeqBAIJ_N(Mat A,Vec xx,Vec zz) 582 { 583 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 584 PetscScalar *x,*z,*xb,*work,*workt,*zarray; 585 MatScalar *v; 586 PetscErrorCode ierr; 587 PetscInt mbs=a->mbs,i,*idx,*ii,bs=A->bs,j,n,bs2=a->bs2; 588 PetscInt ncols,k,*rindex=PETSC_NULL; 589 PetscTruth usecprow=a->compressedrow.use; 590 591 PetscFunctionBegin; 592 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 593 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 594 595 idx = a->j; 596 v = a->a; 597 if (usecprow){ 598 mbs = a->compressedrow.nrows; 599 ii = a->compressedrow.i; 600 rindex = a->compressedrow.rindex; 601 } else { 602 mbs = a->mbs; 603 ii = a->i; 604 z = zarray; 605 } 606 607 if (!a->mult_work) { 608 k = PetscMax(A->m,A->n); 609 ierr = PetscMalloc((k+1)*sizeof(PetscScalar),&a->mult_work);CHKERRQ(ierr); 610 } 611 work = a->mult_work; 612 for (i=0; i<mbs; i++) { 613 n = ii[1] - ii[0]; ii++; 614 ncols = n*bs; 615 workt = work; 616 for (j=0; j<n; j++) { 617 xb = x + bs*(*idx++); 618 for (k=0; k<bs; k++) workt[k] = xb[k]; 619 workt += bs; 620 } 621 if (usecprow) z = zarray + bs*rindex[i]; 622 Kernel_w_gets_Ar_times_v(bs,ncols,work,v,z); 623 /* LAgemv_("N",&bs,&ncols,&_DOne,v,&bs,work,&_One,&_DZero,z,&_One); */ 624 v += n*bs2; 625 if (!usecprow) z += bs; 626 } 627 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 628 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 629 PetscLogFlops(2*a->nz*bs2 - bs*mbs); 630 PetscFunctionReturn(0); 631 } 632 633 #undef __FUNCT__ 634 #define __FUNCT__ "MatMultAdd_SeqBAIJ_1" 635 PetscErrorCode MatMultAdd_SeqBAIJ_1(Mat A,Vec xx,Vec yy,Vec zz) 636 { 637 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 638 PetscScalar *x,*y,*z,sum,*yarray,*zarray; 639 MatScalar *v; 640 PetscErrorCode ierr; 641 PetscInt mbs,i,*idx,*ii,n,*rindex=PETSC_NULL; 642 PetscTruth usecprow=a->compressedrow.use; 643 644 PetscFunctionBegin; 645 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 646 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 647 if (zz != yy) { 648 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 649 } else { 650 zarray = yarray; 651 } 652 653 idx = a->j; 654 v = a->a; 655 if (usecprow){ 656 mbs = a->compressedrow.nrows; 657 ii = a->compressedrow.i; 658 rindex = a->compressedrow.rindex; 659 } else { 660 mbs = a->mbs; 661 ii = a->i; 662 y = yarray; 663 z = zarray; 664 } 665 666 for (i=0; i<mbs; i++) { 667 n = ii[1] - ii[0]; ii++; 668 if (usecprow){ 669 z = zarray + rindex[i]; 670 y = yarray + rindex[i]; 671 } 672 sum = y[0]; 673 while (n--) sum += *v++ * x[*idx++]; 674 z[0] = sum; 675 if (!usecprow){ 676 z++; y++; 677 } 678 } 679 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 680 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 681 if (zz != yy) { 682 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 683 } 684 PetscLogFlops(2*a->nz); 685 PetscFunctionReturn(0); 686 } 687 688 #undef __FUNCT__ 689 #define __FUNCT__ "MatMultAdd_SeqBAIJ_2" 690 PetscErrorCode MatMultAdd_SeqBAIJ_2(Mat A,Vec xx,Vec yy,Vec zz) 691 { 692 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 693 PetscScalar *x,*y,*z,*xb,sum1,sum2; 694 PetscScalar x1,x2,*yarray,*zarray; 695 MatScalar *v; 696 PetscErrorCode ierr; 697 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 698 PetscTruth usecprow=a->compressedrow.use; 699 700 PetscFunctionBegin; 701 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 702 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 703 if (zz != yy) { 704 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 705 } else { 706 zarray = yarray; 707 } 708 709 idx = a->j; 710 v = a->a; 711 if (usecprow){ 712 mbs = a->compressedrow.nrows; 713 ii = a->compressedrow.i; 714 rindex = a->compressedrow.rindex; 715 } else { 716 mbs = a->mbs; 717 ii = a->i; 718 y = yarray; 719 z = zarray; 720 } 721 722 for (i=0; i<mbs; i++) { 723 n = ii[1] - ii[0]; ii++; 724 if (usecprow){ 725 z = zarray + 2*rindex[i]; 726 y = yarray + 2*rindex[i]; 727 } 728 sum1 = y[0]; sum2 = y[1]; 729 for (j=0; j<n; j++) { 730 xb = x + 2*(*idx++); x1 = xb[0]; x2 = xb[1]; 731 sum1 += v[0]*x1 + v[2]*x2; 732 sum2 += v[1]*x1 + v[3]*x2; 733 v += 4; 734 } 735 z[0] = sum1; z[1] = sum2; 736 if (!usecprow){ 737 z += 2; y += 2; 738 } 739 } 740 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 741 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 742 if (zz != yy) { 743 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 744 } 745 PetscLogFlops(4*a->nz); 746 PetscFunctionReturn(0); 747 } 748 749 #undef __FUNCT__ 750 #define __FUNCT__ "MatMultAdd_SeqBAIJ_3" 751 PetscErrorCode MatMultAdd_SeqBAIJ_3(Mat A,Vec xx,Vec yy,Vec zz) 752 { 753 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 754 PetscScalar *x,*y,*z,*xb,sum1,sum2,sum3,x1,x2,x3,*yarray,*zarray; 755 MatScalar *v; 756 PetscErrorCode ierr; 757 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 758 PetscTruth usecprow=a->compressedrow.use; 759 760 PetscFunctionBegin; 761 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 762 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 763 if (zz != yy) { 764 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 765 } else { 766 zarray = yarray; 767 } 768 769 idx = a->j; 770 v = a->a; 771 if (usecprow){ 772 mbs = a->compressedrow.nrows; 773 ii = a->compressedrow.i; 774 rindex = a->compressedrow.rindex; 775 } else { 776 mbs = a->mbs; 777 ii = a->i; 778 y = yarray; 779 z = zarray; 780 } 781 782 for (i=0; i<mbs; i++) { 783 n = ii[1] - ii[0]; ii++; 784 if (usecprow){ 785 z = zarray + 3*rindex[i]; 786 y = yarray + 3*rindex[i]; 787 } 788 sum1 = y[0]; sum2 = y[1]; sum3 = y[2]; 789 for (j=0; j<n; j++) { 790 xb = x + 3*(*idx++); x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 791 sum1 += v[0]*x1 + v[3]*x2 + v[6]*x3; 792 sum2 += v[1]*x1 + v[4]*x2 + v[7]*x3; 793 sum3 += v[2]*x1 + v[5]*x2 + v[8]*x3; 794 v += 9; 795 } 796 z[0] = sum1; z[1] = sum2; z[2] = sum3; 797 if (!usecprow){ 798 z += 3; y += 3; 799 } 800 } 801 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 802 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 803 if (zz != yy) { 804 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 805 } 806 PetscLogFlops(18*a->nz); 807 PetscFunctionReturn(0); 808 } 809 810 #undef __FUNCT__ 811 #define __FUNCT__ "MatMultAdd_SeqBAIJ_4" 812 PetscErrorCode MatMultAdd_SeqBAIJ_4(Mat A,Vec xx,Vec yy,Vec zz) 813 { 814 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 815 PetscScalar *x,*y,*z,*xb,sum1,sum2,sum3,sum4,x1,x2,x3,x4,*yarray,*zarray; 816 MatScalar *v; 817 PetscErrorCode ierr; 818 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 819 PetscTruth usecprow=a->compressedrow.use; 820 821 PetscFunctionBegin; 822 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 823 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 824 if (zz != yy) { 825 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 826 } else { 827 zarray = yarray; 828 } 829 830 idx = a->j; 831 v = a->a; 832 if (usecprow){ 833 mbs = a->compressedrow.nrows; 834 ii = a->compressedrow.i; 835 rindex = a->compressedrow.rindex; 836 } else { 837 mbs = a->mbs; 838 ii = a->i; 839 y = yarray; 840 z = zarray; 841 } 842 843 for (i=0; i<mbs; i++) { 844 n = ii[1] - ii[0]; ii++; 845 if (usecprow){ 846 z = zarray + 4*rindex[i]; 847 y = yarray + 4*rindex[i]; 848 } 849 sum1 = y[0]; sum2 = y[1]; sum3 = y[2]; sum4 = y[3]; 850 for (j=0; j<n; j++) { 851 xb = x + 4*(*idx++); 852 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; 853 sum1 += v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 854 sum2 += v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 855 sum3 += v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 856 sum4 += v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 857 v += 16; 858 } 859 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; 860 if (!usecprow){ 861 z += 4; y += 4; 862 } 863 } 864 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 865 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 866 if (zz != yy) { 867 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 868 } 869 PetscLogFlops(32*a->nz); 870 PetscFunctionReturn(0); 871 } 872 873 #undef __FUNCT__ 874 #define __FUNCT__ "MatMultAdd_SeqBAIJ_5" 875 PetscErrorCode MatMultAdd_SeqBAIJ_5(Mat A,Vec xx,Vec yy,Vec zz) 876 { 877 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 878 PetscScalar *x,*y,*z,*xb,sum1,sum2,sum3,sum4,sum5,x1,x2,x3,x4,x5; 879 PetscScalar *yarray,*zarray; 880 MatScalar *v; 881 PetscErrorCode ierr; 882 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 883 PetscTruth usecprow=a->compressedrow.use; 884 885 PetscFunctionBegin; 886 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 887 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 888 if (zz != yy) { 889 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 890 } else { 891 zarray = yarray; 892 } 893 894 idx = a->j; 895 v = a->a; 896 if (usecprow){ 897 mbs = a->compressedrow.nrows; 898 ii = a->compressedrow.i; 899 rindex = a->compressedrow.rindex; 900 } else { 901 mbs = a->mbs; 902 ii = a->i; 903 y = yarray; 904 z = zarray; 905 } 906 907 for (i=0; i<mbs; i++) { 908 n = ii[1] - ii[0]; ii++; 909 if (usecprow){ 910 z = zarray + 5*rindex[i]; 911 y = yarray + 5*rindex[i]; 912 } 913 sum1 = y[0]; sum2 = y[1]; sum3 = y[2]; sum4 = y[3]; sum5 = y[4]; 914 for (j=0; j<n; j++) { 915 xb = x + 5*(*idx++); 916 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; 917 sum1 += v[0]*x1 + v[5]*x2 + v[10]*x3 + v[15]*x4 + v[20]*x5; 918 sum2 += v[1]*x1 + v[6]*x2 + v[11]*x3 + v[16]*x4 + v[21]*x5; 919 sum3 += v[2]*x1 + v[7]*x2 + v[12]*x3 + v[17]*x4 + v[22]*x5; 920 sum4 += v[3]*x1 + v[8]*x2 + v[13]*x3 + v[18]*x4 + v[23]*x5; 921 sum5 += v[4]*x1 + v[9]*x2 + v[14]*x3 + v[19]*x4 + v[24]*x5; 922 v += 25; 923 } 924 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; 925 if (!usecprow){ 926 z += 5; y += 5; 927 } 928 } 929 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 930 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 931 if (zz != yy) { 932 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 933 } 934 PetscLogFlops(50*a->nz); 935 PetscFunctionReturn(0); 936 } 937 #undef __FUNCT__ 938 #define __FUNCT__ "MatMultAdd_SeqBAIJ_6" 939 PetscErrorCode MatMultAdd_SeqBAIJ_6(Mat A,Vec xx,Vec yy,Vec zz) 940 { 941 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 942 PetscScalar *x,*y,*z,*xb,sum1,sum2,sum3,sum4,sum5,sum6; 943 PetscScalar x1,x2,x3,x4,x5,x6,*yarray,*zarray; 944 MatScalar *v; 945 PetscErrorCode ierr; 946 PetscInt mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 947 PetscTruth usecprow=a->compressedrow.use; 948 949 PetscFunctionBegin; 950 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 951 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 952 if (zz != yy) { 953 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 954 } else { 955 zarray = yarray; 956 } 957 958 idx = a->j; 959 v = a->a; 960 if (usecprow){ 961 mbs = a->compressedrow.nrows; 962 ii = a->compressedrow.i; 963 rindex = a->compressedrow.rindex; 964 } else { 965 mbs = a->mbs; 966 ii = a->i; 967 y = yarray; 968 z = zarray; 969 } 970 971 for (i=0; i<mbs; i++) { 972 n = ii[1] - ii[0]; ii++; 973 if (usecprow){ 974 z = zarray + 6*rindex[i]; 975 y = yarray + 6*rindex[i]; 976 } 977 sum1 = y[0]; sum2 = y[1]; sum3 = y[2]; sum4 = y[3]; sum5 = y[4]; sum6 = y[5]; 978 for (j=0; j<n; j++) { 979 xb = x + 6*(*idx++); 980 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; 981 sum1 += v[0]*x1 + v[6]*x2 + v[12]*x3 + v[18]*x4 + v[24]*x5 + v[30]*x6; 982 sum2 += v[1]*x1 + v[7]*x2 + v[13]*x3 + v[19]*x4 + v[25]*x5 + v[31]*x6; 983 sum3 += v[2]*x1 + v[8]*x2 + v[14]*x3 + v[20]*x4 + v[26]*x5 + v[32]*x6; 984 sum4 += v[3]*x1 + v[9]*x2 + v[15]*x3 + v[21]*x4 + v[27]*x5 + v[33]*x6; 985 sum5 += v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4 + v[28]*x5 + v[34]*x6; 986 sum6 += v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4 + v[29]*x5 + v[35]*x6; 987 v += 36; 988 } 989 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; z[5] = sum6; 990 if (!usecprow){ 991 z += 6; y += 6; 992 } 993 } 994 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 995 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 996 if (zz != yy) { 997 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 998 } 999 PetscLogFlops(72*a->nz); 1000 PetscFunctionReturn(0); 1001 } 1002 1003 #undef __FUNCT__ 1004 #define __FUNCT__ "MatMultAdd_SeqBAIJ_7" 1005 PetscErrorCode MatMultAdd_SeqBAIJ_7(Mat A,Vec xx,Vec yy,Vec zz) 1006 { 1007 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1008 PetscScalar *x,*y,*z,*xb,sum1,sum2,sum3,sum4,sum5,sum6,sum7; 1009 PetscScalar x1,x2,x3,x4,x5,x6,x7,*yarray,*zarray; 1010 MatScalar *v; 1011 PetscErrorCode ierr; 1012 PetscInt mbs=a->mbs,i,*idx,*ii,j,n,*rindex=PETSC_NULL; 1013 PetscTruth usecprow=a->compressedrow.use; 1014 1015 PetscFunctionBegin; 1016 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1017 ierr = VecGetArray(yy,&yarray);CHKERRQ(ierr); 1018 if (zz != yy) { 1019 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 1020 } else { 1021 zarray = yarray; 1022 } 1023 1024 idx = a->j; 1025 v = a->a; 1026 if (usecprow){ 1027 mbs = a->compressedrow.nrows; 1028 ii = a->compressedrow.i; 1029 rindex = a->compressedrow.rindex; 1030 } else { 1031 mbs = a->mbs; 1032 ii = a->i; 1033 y = yarray; 1034 z = zarray; 1035 } 1036 1037 for (i=0; i<mbs; i++) { 1038 n = ii[1] - ii[0]; ii++; 1039 if (usecprow){ 1040 z = zarray + 7*rindex[i]; 1041 y = yarray + 7*rindex[i]; 1042 } 1043 sum1 = y[0]; sum2 = y[1]; sum3 = y[2]; sum4 = y[3]; sum5 = y[4]; sum6 = y[5]; sum7 = y[6]; 1044 for (j=0; j<n; j++) { 1045 xb = x + 7*(*idx++); 1046 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; x7 = xb[6]; 1047 sum1 += v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 1048 sum2 += v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 1049 sum3 += v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 1050 sum4 += v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 1051 sum5 += v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 1052 sum6 += v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 1053 sum7 += v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 1054 v += 49; 1055 } 1056 z[0] = sum1; z[1] = sum2; z[2] = sum3; z[3] = sum4; z[4] = sum5; z[5] = sum6; z[6] = sum7; 1057 if (!usecprow){ 1058 z += 7; y += 7; 1059 } 1060 } 1061 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1062 ierr = VecRestoreArray(yy,&yarray);CHKERRQ(ierr); 1063 if (zz != yy) { 1064 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 1065 } 1066 PetscLogFlops(98*a->nz); 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "MatMultAdd_SeqBAIJ_N" 1072 PetscErrorCode MatMultAdd_SeqBAIJ_N(Mat A,Vec xx,Vec yy,Vec zz) 1073 { 1074 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1075 PetscScalar *x,*z,*xb,*work,*workt,*y,*zarray; 1076 MatScalar *v; 1077 PetscErrorCode ierr; 1078 PetscInt mbs,i,*idx,*ii,bs=A->bs,j,n,bs2=a->bs2; 1079 PetscInt ncols,k,*rindex=PETSC_NULL; 1080 PetscTruth usecprow=a->compressedrow.use; 1081 1082 PetscFunctionBegin; 1083 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1084 ierr = VecGetArray(zz,&zarray);CHKERRQ(ierr); 1085 if (zz != yy) { 1086 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1087 ierr = PetscMemcpy(zarray,y,yy->n*sizeof(PetscScalar));CHKERRQ(ierr); 1088 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1089 } 1090 1091 idx = a->j; 1092 v = a->a; 1093 if (usecprow){ 1094 mbs = a->compressedrow.nrows; 1095 ii = a->compressedrow.i; 1096 rindex = a->compressedrow.rindex; 1097 } else { 1098 mbs = a->mbs; 1099 ii = a->i; 1100 z = zarray; 1101 } 1102 1103 if (!a->mult_work) { 1104 k = PetscMax(A->m,A->n); 1105 ierr = PetscMalloc((k+1)*sizeof(PetscScalar),&a->mult_work);CHKERRQ(ierr); 1106 } 1107 work = a->mult_work; 1108 for (i=0; i<mbs; i++) { 1109 n = ii[1] - ii[0]; ii++; 1110 ncols = n*bs; 1111 workt = work; 1112 for (j=0; j<n; j++) { 1113 xb = x + bs*(*idx++); 1114 for (k=0; k<bs; k++) workt[k] = xb[k]; 1115 workt += bs; 1116 } 1117 if (usecprow) z = zarray + bs*rindex[i]; 1118 Kernel_w_gets_w_plus_Ar_times_v(bs,ncols,work,v,z); 1119 /* LAgemv_("N",&bs,&ncols,&_DOne,v,&bs,work,&_One,&_DOne,z,&_One); */ 1120 v += n*bs2; 1121 if (!usecprow){ 1122 z += bs; 1123 } 1124 } 1125 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1126 ierr = VecRestoreArray(zz,&zarray);CHKERRQ(ierr); 1127 PetscLogFlops(2*a->nz*bs2); 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "MatMultTranspose_SeqBAIJ" 1133 PetscErrorCode MatMultTranspose_SeqBAIJ(Mat A,Vec xx,Vec zz) 1134 { 1135 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1136 PetscScalar *xg,*zg,*zb,zero = 0.0; 1137 PetscScalar *x,*z,*xb,x1,x2,x3,x4,x5,x6,x7; 1138 MatScalar *v; 1139 PetscInt mbs=a->mbs,i,*idx,*ii,*ai=a->i,rval; 1140 PetscErrorCode ierr; 1141 PetscInt bs=A->bs,j,n,bs2=a->bs2,*ib; 1142 1143 PetscFunctionBegin; 1144 ierr = VecSet(&zero,zz);CHKERRQ(ierr); 1145 ierr = VecGetArray(xx,&xg);CHKERRQ(ierr); x = xg; 1146 ierr = VecGetArray(zz,&zg);CHKERRQ(ierr); z = zg; 1147 1148 idx = a->j; 1149 v = a->a; 1150 ii = a->i; 1151 xb = x; 1152 switch (bs) { 1153 case 1: 1154 for (i=0; i<mbs; i++) { 1155 n = ii[1] - ii[0]; ii++; 1156 x1 = xb[0]; 1157 ib = idx + ai[i]; 1158 for (j=0; j<n; j++) { 1159 rval = ib[j]; 1160 z[rval] += *v * x1; 1161 v++; 1162 } 1163 xb++; 1164 } 1165 break; 1166 case 2: 1167 for (i=0; i<mbs; i++) { 1168 n = ii[1] - ii[0]; ii++; 1169 x1 = xb[0]; x2 = xb[1]; 1170 ib = idx + ai[i]; 1171 for (j=0; j<n; j++) { 1172 rval = ib[j]*2; 1173 z[rval++] += v[0]*x1 + v[1]*x2; 1174 z[rval] += v[2]*x1 + v[3]*x2; 1175 v += 4; 1176 } 1177 xb += 2; 1178 } 1179 break; 1180 case 3: 1181 for (i=0; i<mbs; i++) { 1182 n = ii[1] - ii[0]; ii++; 1183 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1184 ib = idx + ai[i]; 1185 for (j=0; j<n; j++) { 1186 rval = ib[j]*3; 1187 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3; 1188 z[rval++] += v[3]*x1 + v[4]*x2 + v[5]*x3; 1189 z[rval] += v[6]*x1 + v[7]*x2 + v[8]*x3; 1190 v += 9; 1191 } 1192 xb += 3; 1193 } 1194 break; 1195 case 4: 1196 for (i=0; i<mbs; i++) { 1197 n = ii[1] - ii[0]; ii++; 1198 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; 1199 ib = idx + ai[i]; 1200 for (j=0; j<n; j++) { 1201 rval = ib[j]*4; 1202 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4; 1203 z[rval++] += v[4]*x1 + v[5]*x2 + v[6]*x3 + v[7]*x4; 1204 z[rval++] += v[8]*x1 + v[9]*x2 + v[10]*x3 + v[11]*x4; 1205 z[rval] += v[12]*x1 + v[13]*x2 + v[14]*x3 + v[15]*x4; 1206 v += 16; 1207 } 1208 xb += 4; 1209 } 1210 break; 1211 case 5: 1212 for (i=0; i<mbs; i++) { 1213 n = ii[1] - ii[0]; ii++; 1214 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1215 x4 = xb[3]; x5 = xb[4]; 1216 ib = idx + ai[i]; 1217 for (j=0; j<n; j++) { 1218 rval = ib[j]*5; 1219 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5; 1220 z[rval++] += v[5]*x1 + v[6]*x2 + v[7]*x3 + v[8]*x4 + v[9]*x5; 1221 z[rval++] += v[10]*x1 + v[11]*x2 + v[12]*x3 + v[13]*x4 + v[14]*x5; 1222 z[rval++] += v[15]*x1 + v[16]*x2 + v[17]*x3 + v[18]*x4 + v[19]*x5; 1223 z[rval] += v[20]*x1 + v[21]*x2 + v[22]*x3 + v[23]*x4 + v[24]*x5; 1224 v += 25; 1225 } 1226 xb += 5; 1227 } 1228 break; 1229 case 6: 1230 for (i=0; i<mbs; i++) { 1231 n = ii[1] - ii[0]; ii++; 1232 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1233 x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; 1234 ib = idx + ai[i]; 1235 for (j=0; j<n; j++) { 1236 rval = ib[j]*6; 1237 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5 + v[5]*x6; 1238 z[rval++] += v[6]*x1 + v[7]*x2 + v[8]*x3 + v[9]*x4 + v[10]*x5 + v[11]*x6; 1239 z[rval++] += v[12]*x1 + v[13]*x2 + v[14]*x3 + v[15]*x4 + v[16]*x5 + v[17]*x6; 1240 z[rval++] += v[18]*x1 + v[19]*x2 + v[20]*x3 + v[21]*x4 + v[22]*x5 + v[23]*x6; 1241 z[rval++] += v[24]*x1 + v[25]*x2 + v[26]*x3 + v[27]*x4 + v[28]*x5 + v[29]*x6; 1242 z[rval] += v[30]*x1 + v[31]*x2 + v[32]*x3 + v[33]*x4 + v[34]*x5 + v[35]*x6; 1243 v += 36; 1244 } 1245 xb += 6; 1246 } 1247 break; 1248 case 7: 1249 for (i=0; i<mbs; i++) { 1250 n = ii[1] - ii[0]; ii++; 1251 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1252 x4 = xb[3]; x5 = xb[4]; x6 = xb[5]; x7 = xb[6]; 1253 ib = idx + ai[i]; 1254 for (j=0; j<n; j++) { 1255 rval = ib[j]*7; 1256 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5 + v[5]*x6 + v[6]*x7; 1257 z[rval++] += v[7]*x1 + v[8]*x2 + v[9]*x3 + v[10]*x4 + v[11]*x5 + v[12]*x6 + v[13]*x7; 1258 z[rval++] += v[14]*x1 + v[15]*x2 + v[16]*x3 + v[17]*x4 + v[18]*x5 + v[19]*x6 + v[20]*x7; 1259 z[rval++] += v[21]*x1 + v[22]*x2 + v[23]*x3 + v[24]*x4 + v[25]*x5 + v[26]*x6 + v[27]*x7; 1260 z[rval++] += v[28]*x1 + v[29]*x2 + v[30]*x3 + v[31]*x4 + v[32]*x5 + v[33]*x6 + v[34]*x7; 1261 z[rval++] += v[35]*x1 + v[36]*x2 + v[37]*x3 + v[38]*x4 + v[39]*x5 + v[40]*x6 + v[41]*x7; 1262 z[rval] += v[42]*x1 + v[43]*x2 + v[44]*x3 + v[45]*x4 + v[46]*x5 + v[47]*x6 + v[48]*x7; 1263 v += 49; 1264 } 1265 xb += 7; 1266 } 1267 break; 1268 default: { /* block sizes larger then 7 by 7 are handled by BLAS */ 1269 PetscInt ncols,k; 1270 PetscScalar *work,*workt; 1271 1272 if (!a->mult_work) { 1273 k = PetscMax(A->m,A->n); 1274 ierr = PetscMalloc((k+1)*sizeof(PetscScalar),&a->mult_work);CHKERRQ(ierr); 1275 } 1276 work = a->mult_work; 1277 for (i=0; i<mbs; i++) { 1278 n = ii[1] - ii[0]; ii++; 1279 ncols = n*bs; 1280 ierr = PetscMemzero(work,ncols*sizeof(PetscScalar));CHKERRQ(ierr); 1281 Kernel_w_gets_w_plus_trans_Ar_times_v(bs,ncols,x,v,work); 1282 /* LAgemv_("T",&bs,&ncols,&_DOne,v,&bs,x,&_One,&_DOne,work,&_One); */ 1283 v += n*bs2; 1284 x += bs; 1285 workt = work; 1286 for (j=0; j<n; j++) { 1287 zb = z + bs*(*idx++); 1288 for (k=0; k<bs; k++) zb[k] += workt[k] ; 1289 workt += bs; 1290 } 1291 } 1292 } 1293 } 1294 ierr = VecRestoreArray(xx,&xg);CHKERRQ(ierr); 1295 ierr = VecRestoreArray(zz,&zg);CHKERRQ(ierr); 1296 PetscLogFlops(2*a->nz*a->bs2 - A->n); 1297 PetscFunctionReturn(0); 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "MatMultTransposeAdd_SeqBAIJ" 1302 PetscErrorCode MatMultTransposeAdd_SeqBAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1303 1304 { 1305 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1306 PetscScalar *xg,*zg,*zb,*x,*z,*xb,x1,x2,x3,x4,x5; 1307 MatScalar *v; 1308 PetscErrorCode ierr; 1309 PetscInt mbs=a->mbs,i,*idx,*ii,*ai=a->i,rval,bs=A->bs,j,n,bs2=a->bs2,*ib; 1310 1311 PetscFunctionBegin; 1312 if (yy != zz) { ierr = VecCopy(yy,zz);CHKERRQ(ierr); } 1313 ierr = VecGetArray(xx,&xg);CHKERRQ(ierr); x = xg; 1314 ierr = VecGetArray(zz,&zg);CHKERRQ(ierr); z = zg; 1315 1316 1317 idx = a->j; 1318 v = a->a; 1319 ii = a->i; 1320 xb = x; 1321 1322 switch (bs) { 1323 case 1: 1324 for (i=0; i<mbs; i++) { 1325 n = ii[1] - ii[0]; ii++; 1326 x1 = xb[0]; 1327 ib = idx + ai[i]; 1328 for (j=0; j<n; j++) { 1329 rval = ib[j]; 1330 z[rval] += *v * x1; 1331 v++; 1332 } 1333 xb++; 1334 } 1335 break; 1336 case 2: 1337 for (i=0; i<mbs; i++) { 1338 n = ii[1] - ii[0]; ii++; 1339 x1 = xb[0]; x2 = xb[1]; 1340 ib = idx + ai[i]; 1341 for (j=0; j<n; j++) { 1342 rval = ib[j]*2; 1343 z[rval++] += v[0]*x1 + v[1]*x2; 1344 z[rval++] += v[2]*x1 + v[3]*x2; 1345 v += 4; 1346 } 1347 xb += 2; 1348 } 1349 break; 1350 case 3: 1351 for (i=0; i<mbs; i++) { 1352 n = ii[1] - ii[0]; ii++; 1353 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1354 ib = idx + ai[i]; 1355 for (j=0; j<n; j++) { 1356 rval = ib[j]*3; 1357 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3; 1358 z[rval++] += v[3]*x1 + v[4]*x2 + v[5]*x3; 1359 z[rval++] += v[6]*x1 + v[7]*x2 + v[8]*x3; 1360 v += 9; 1361 } 1362 xb += 3; 1363 } 1364 break; 1365 case 4: 1366 for (i=0; i<mbs; i++) { 1367 n = ii[1] - ii[0]; ii++; 1368 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; 1369 ib = idx + ai[i]; 1370 for (j=0; j<n; j++) { 1371 rval = ib[j]*4; 1372 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4; 1373 z[rval++] += v[4]*x1 + v[5]*x2 + v[6]*x3 + v[7]*x4; 1374 z[rval++] += v[8]*x1 + v[9]*x2 + v[10]*x3 + v[11]*x4; 1375 z[rval++] += v[12]*x1 + v[13]*x2 + v[14]*x3 + v[15]*x4; 1376 v += 16; 1377 } 1378 xb += 4; 1379 } 1380 break; 1381 case 5: 1382 for (i=0; i<mbs; i++) { 1383 n = ii[1] - ii[0]; ii++; 1384 x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; 1385 x4 = xb[3]; x5 = xb[4]; 1386 ib = idx + ai[i]; 1387 for (j=0; j<n; j++) { 1388 rval = ib[j]*5; 1389 z[rval++] += v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5; 1390 z[rval++] += v[5]*x1 + v[6]*x2 + v[7]*x3 + v[8]*x4 + v[9]*x5; 1391 z[rval++] += v[10]*x1 + v[11]*x2 + v[12]*x3 + v[13]*x4 + v[14]*x5; 1392 z[rval++] += v[15]*x1 + v[16]*x2 + v[17]*x3 + v[18]*x4 + v[19]*x5; 1393 z[rval++] += v[20]*x1 + v[21]*x2 + v[22]*x3 + v[23]*x4 + v[24]*x5; 1394 v += 25; 1395 } 1396 xb += 5; 1397 } 1398 break; 1399 default: { /* block sizes larger then 5 by 5 are handled by BLAS */ 1400 PetscInt ncols,k; 1401 PetscScalar *work,*workt; 1402 1403 if (!a->mult_work) { 1404 k = PetscMax(A->m,A->n); 1405 ierr = PetscMalloc((k+1)*sizeof(PetscScalar),&a->mult_work);CHKERRQ(ierr); 1406 } 1407 work = a->mult_work; 1408 for (i=0; i<mbs; i++) { 1409 n = ii[1] - ii[0]; ii++; 1410 ncols = n*bs; 1411 ierr = PetscMemzero(work,ncols*sizeof(PetscScalar));CHKERRQ(ierr); 1412 Kernel_w_gets_w_plus_trans_Ar_times_v(bs,ncols,x,v,work); 1413 /* LAgemv_("T",&bs,&ncols,&_DOne,v,&bs,x,&_One,&_DOne,work,&_One); */ 1414 v += n*bs2; 1415 x += bs; 1416 workt = work; 1417 for (j=0; j<n; j++) { 1418 zb = z + bs*(*idx++); 1419 for (k=0; k<bs; k++) zb[k] += workt[k] ; 1420 workt += bs; 1421 } 1422 } 1423 } 1424 } 1425 ierr = VecRestoreArray(xx,&xg);CHKERRQ(ierr); 1426 ierr = VecRestoreArray(zz,&zg);CHKERRQ(ierr); 1427 PetscLogFlops(2*a->nz*a->bs2); 1428 PetscFunctionReturn(0); 1429 } 1430 1431 #undef __FUNCT__ 1432 #define __FUNCT__ "MatScale_SeqBAIJ" 1433 PetscErrorCode MatScale_SeqBAIJ(const PetscScalar *alpha,Mat inA) 1434 { 1435 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)inA->data; 1436 PetscInt totalnz = a->bs2*a->nz; 1437 #if defined(PETSC_USE_MAT_SINGLE) 1438 PetscInt i; 1439 #else 1440 PetscBLASInt tnz = (PetscBLASInt) totalnz,one = 1; 1441 #endif 1442 1443 PetscFunctionBegin; 1444 #if defined(PETSC_USE_MAT_SINGLE) 1445 for (i=0; i<totalnz; i++) a->a[i] *= *alpha; 1446 #else 1447 BLscal_(&tnz,(PetscScalar*)alpha,a->a,&one); 1448 #endif 1449 PetscLogFlops(totalnz); 1450 PetscFunctionReturn(0); 1451 } 1452 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "MatNorm_SeqBAIJ" 1455 PetscErrorCode MatNorm_SeqBAIJ(Mat A,NormType type,PetscReal *norm) 1456 { 1457 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1458 MatScalar *v = a->a; 1459 PetscReal sum = 0.0; 1460 PetscInt i,j,k,bs = A->bs,nz=a->nz,bs2=a->bs2,k1; 1461 1462 PetscFunctionBegin; 1463 if (type == NORM_FROBENIUS) { 1464 for (i=0; i< bs2*nz; i++) { 1465 #if defined(PETSC_USE_COMPLEX) 1466 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1467 #else 1468 sum += (*v)*(*v); v++; 1469 #endif 1470 } 1471 *norm = sqrt(sum); 1472 } else if (type == NORM_INFINITY) { /* maximum row sum */ 1473 *norm = 0.0; 1474 for (k=0; k<bs; k++) { 1475 for (j=0; j<a->mbs; j++) { 1476 v = a->a + bs2*a->i[j] + k; 1477 sum = 0.0; 1478 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1479 for (k1=0; k1<bs; k1++){ 1480 sum += PetscAbsScalar(*v); 1481 v += bs; 1482 } 1483 } 1484 if (sum > *norm) *norm = sum; 1485 } 1486 } 1487 } else { 1488 SETERRQ(PETSC_ERR_SUP,"No support for this norm yet"); 1489 } 1490 PetscFunctionReturn(0); 1491 } 1492 1493 1494 #undef __FUNCT__ 1495 #define __FUNCT__ "MatEqual_SeqBAIJ" 1496 PetscErrorCode MatEqual_SeqBAIJ(Mat A,Mat B,PetscTruth* flg) 1497 { 1498 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data,*b = (Mat_SeqBAIJ *)B->data; 1499 PetscErrorCode ierr; 1500 1501 PetscFunctionBegin; 1502 /* If the matrix/block dimensions are not equal, or no of nonzeros or shift */ 1503 if ((A->m != B->m) || (A->n != B->n) || (A->bs != B->bs)|| (a->nz != b->nz)) { 1504 *flg = PETSC_FALSE; 1505 PetscFunctionReturn(0); 1506 } 1507 1508 /* if the a->i are the same */ 1509 ierr = PetscMemcmp(a->i,b->i,(a->mbs+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 1510 if (*flg == PETSC_FALSE) { 1511 PetscFunctionReturn(0); 1512 } 1513 1514 /* if a->j are the same */ 1515 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 1516 if (*flg == PETSC_FALSE) { 1517 PetscFunctionReturn(0); 1518 } 1519 /* if a->a are the same */ 1520 ierr = PetscMemcmp(a->a,b->a,(a->nz)*(A->bs)*(B->bs)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 1521 PetscFunctionReturn(0); 1522 1523 } 1524 1525 #undef __FUNCT__ 1526 #define __FUNCT__ "MatGetDiagonal_SeqBAIJ" 1527 PetscErrorCode MatGetDiagonal_SeqBAIJ(Mat A,Vec v) 1528 { 1529 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1530 PetscErrorCode ierr; 1531 PetscInt i,j,k,n,row,bs,*ai,*aj,ambs,bs2; 1532 PetscScalar *x,zero = 0.0; 1533 MatScalar *aa,*aa_j; 1534 1535 PetscFunctionBegin; 1536 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1537 bs = A->bs; 1538 aa = a->a; 1539 ai = a->i; 1540 aj = a->j; 1541 ambs = a->mbs; 1542 bs2 = a->bs2; 1543 1544 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1545 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1546 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1547 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1548 for (i=0; i<ambs; i++) { 1549 for (j=ai[i]; j<ai[i+1]; j++) { 1550 if (aj[j] == i) { 1551 row = i*bs; 1552 aa_j = aa+j*bs2; 1553 for (k=0; k<bs2; k+=(bs+1),row++) x[row] = aa_j[k]; 1554 break; 1555 } 1556 } 1557 } 1558 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1559 PetscFunctionReturn(0); 1560 } 1561 1562 #undef __FUNCT__ 1563 #define __FUNCT__ "MatDiagonalScale_SeqBAIJ" 1564 PetscErrorCode MatDiagonalScale_SeqBAIJ(Mat A,Vec ll,Vec rr) 1565 { 1566 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1567 PetscScalar *l,*r,x,*li,*ri; 1568 MatScalar *aa,*v; 1569 PetscErrorCode ierr; 1570 PetscInt i,j,k,lm,rn,M,m,n,*ai,*aj,mbs,tmp,bs,bs2; 1571 1572 PetscFunctionBegin; 1573 ai = a->i; 1574 aj = a->j; 1575 aa = a->a; 1576 m = A->m; 1577 n = A->n; 1578 bs = A->bs; 1579 mbs = a->mbs; 1580 bs2 = a->bs2; 1581 if (ll) { 1582 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1583 ierr = VecGetLocalSize(ll,&lm);CHKERRQ(ierr); 1584 if (lm != m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1585 for (i=0; i<mbs; i++) { /* for each block row */ 1586 M = ai[i+1] - ai[i]; 1587 li = l + i*bs; 1588 v = aa + bs2*ai[i]; 1589 for (j=0; j<M; j++) { /* for each block */ 1590 for (k=0; k<bs2; k++) { 1591 (*v++) *= li[k%bs]; 1592 } 1593 } 1594 } 1595 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1596 PetscLogFlops(a->nz); 1597 } 1598 1599 if (rr) { 1600 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1601 ierr = VecGetLocalSize(rr,&rn);CHKERRQ(ierr); 1602 if (rn != n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1603 for (i=0; i<mbs; i++) { /* for each block row */ 1604 M = ai[i+1] - ai[i]; 1605 v = aa + bs2*ai[i]; 1606 for (j=0; j<M; j++) { /* for each block */ 1607 ri = r + bs*aj[ai[i]+j]; 1608 for (k=0; k<bs; k++) { 1609 x = ri[k]; 1610 for (tmp=0; tmp<bs; tmp++) (*v++) *= x; 1611 } 1612 } 1613 } 1614 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1615 PetscLogFlops(a->nz); 1616 } 1617 PetscFunctionReturn(0); 1618 } 1619 1620 1621 #undef __FUNCT__ 1622 #define __FUNCT__ "MatGetInfo_SeqBAIJ" 1623 PetscErrorCode MatGetInfo_SeqBAIJ(Mat A,MatInfoType flag,MatInfo *info) 1624 { 1625 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1626 1627 PetscFunctionBegin; 1628 info->rows_global = (double)A->m; 1629 info->columns_global = (double)A->n; 1630 info->rows_local = (double)A->m; 1631 info->columns_local = (double)A->n; 1632 info->block_size = a->bs2; 1633 info->nz_allocated = a->maxnz; 1634 info->nz_used = a->bs2*a->nz; 1635 info->nz_unneeded = (double)(info->nz_allocated - info->nz_used); 1636 info->assemblies = A->num_ass; 1637 info->mallocs = a->reallocs; 1638 info->memory = A->mem; 1639 if (A->factor) { 1640 info->fill_ratio_given = A->info.fill_ratio_given; 1641 info->fill_ratio_needed = A->info.fill_ratio_needed; 1642 info->factor_mallocs = A->info.factor_mallocs; 1643 } else { 1644 info->fill_ratio_given = 0; 1645 info->fill_ratio_needed = 0; 1646 info->factor_mallocs = 0; 1647 } 1648 PetscFunctionReturn(0); 1649 } 1650 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatZeroEntries_SeqBAIJ" 1654 PetscErrorCode MatZeroEntries_SeqBAIJ(Mat A) 1655 { 1656 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1657 PetscErrorCode ierr; 1658 1659 PetscFunctionBegin; 1660 ierr = PetscMemzero(a->a,a->bs2*a->i[a->mbs]*sizeof(MatScalar));CHKERRQ(ierr); 1661 PetscFunctionReturn(0); 1662 } 1663