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