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