1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 8 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 9 #include <petscblaslapack.h> 10 #include <petscbt.h> 11 #include <petsc-private/kernels/blocktranspose.h> 12 #if defined(PETSC_THREADCOMM_ACTIVE) 13 #include <petscthreadcomm.h> 14 #endif 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 18 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 19 { 20 PetscErrorCode ierr; 21 PetscInt i,m,n; 22 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 23 24 PetscFunctionBegin; 25 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 26 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 27 if (type == NORM_2) { 28 for (i=0; i<aij->i[m]; i++) { 29 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 30 } 31 } else if (type == NORM_1) { 32 for (i=0; i<aij->i[m]; i++) { 33 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 34 } 35 } else if (type == NORM_INFINITY) { 36 for (i=0; i<aij->i[m]; i++) { 37 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 38 } 39 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 40 41 if (type == NORM_2) { 42 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatFindOffBlockDiagonalEntries_SeqAIJ" 49 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 50 { 51 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 52 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 53 const PetscInt *jj = a->j,*ii = a->i; 54 PetscInt *rows; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 for (i=0; i<m; i++) { 59 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 60 cnt++; 61 } 62 } 63 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 64 cnt = 0; 65 for (i=0; i<m; i++) { 66 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 67 rows[cnt] = i; 68 cnt++; 69 } 70 } 71 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ_Private" 77 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 78 { 79 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 80 const MatScalar *aa = a->a; 81 PetscInt i,m=A->rmap->n,cnt = 0; 82 const PetscInt *jj = a->j,*diag; 83 PetscInt *rows; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 88 diag = a->diag; 89 for (i=0; i<m; i++) { 90 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 91 cnt++; 92 } 93 } 94 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 95 cnt = 0; 96 for (i=0; i<m; i++) { 97 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 98 rows[cnt++] = i; 99 } 100 } 101 *nrows = cnt; 102 *zrows = rows; 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ" 108 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 109 { 110 PetscInt nrows,*rows; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 *zrows = NULL; 115 ierr = MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);CHKERRQ(ierr); 116 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatFindNonzeroRows_SeqAIJ" 122 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 123 { 124 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 125 const MatScalar *aa; 126 PetscInt m=A->rmap->n,cnt = 0; 127 const PetscInt *ii; 128 PetscInt n,i,j,*rows; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 *keptrows = 0; 133 ii = a->i; 134 for (i=0; i<m; i++) { 135 n = ii[i+1] - ii[i]; 136 if (!n) { 137 cnt++; 138 goto ok1; 139 } 140 aa = a->a + ii[i]; 141 for (j=0; j<n; j++) { 142 if (aa[j] != 0.0) goto ok1; 143 } 144 cnt++; 145 ok1:; 146 } 147 if (!cnt) PetscFunctionReturn(0); 148 ierr = PetscMalloc1((A->rmap->n-cnt),&rows);CHKERRQ(ierr); 149 cnt = 0; 150 for (i=0; i<m; i++) { 151 n = ii[i+1] - ii[i]; 152 if (!n) continue; 153 aa = a->a + ii[i]; 154 for (j=0; j<n; j++) { 155 if (aa[j] != 0.0) { 156 rows[cnt++] = i; 157 break; 158 } 159 } 160 } 161 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 167 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 168 { 169 PetscErrorCode ierr; 170 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 171 PetscInt i,*diag, m = Y->rmap->n; 172 MatScalar *aa = aij->a; 173 PetscScalar *v; 174 PetscBool missing; 175 176 PetscFunctionBegin; 177 if (Y->assembled) { 178 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);CHKERRQ(ierr); 179 if (!missing) { 180 diag = aij->diag; 181 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 182 if (is == INSERT_VALUES) { 183 for (i=0; i<m; i++) { 184 aa[diag[i]] = v[i]; 185 } 186 } else { 187 for (i=0; i<m; i++) { 188 aa[diag[i]] += v[i]; 189 } 190 } 191 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 195 } 196 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 202 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 203 { 204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 205 PetscErrorCode ierr; 206 PetscInt i,ishift; 207 208 PetscFunctionBegin; 209 *m = A->rmap->n; 210 if (!ia) PetscFunctionReturn(0); 211 ishift = 0; 212 if (symmetric && !A->structurally_symmetric) { 213 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 214 } else if (oshift == 1) { 215 PetscInt *tia; 216 PetscInt nz = a->i[A->rmap->n]; 217 /* malloc space and add 1 to i and j indices */ 218 ierr = PetscMalloc1((A->rmap->n+1),&tia);CHKERRQ(ierr); 219 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 220 *ia = tia; 221 if (ja) { 222 PetscInt *tja; 223 ierr = PetscMalloc1((nz+1),&tja);CHKERRQ(ierr); 224 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 225 *ja = tja; 226 } 227 } else { 228 *ia = a->i; 229 if (ja) *ja = a->j; 230 } 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 236 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 237 { 238 PetscErrorCode ierr; 239 240 PetscFunctionBegin; 241 if (!ia) PetscFunctionReturn(0); 242 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 243 ierr = PetscFree(*ia);CHKERRQ(ierr); 244 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 245 } 246 PetscFunctionReturn(0); 247 } 248 249 #undef __FUNCT__ 250 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 251 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 252 { 253 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 254 PetscErrorCode ierr; 255 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 256 PetscInt nz = a->i[m],row,*jj,mr,col; 257 258 PetscFunctionBegin; 259 *nn = n; 260 if (!ia) PetscFunctionReturn(0); 261 if (symmetric) { 262 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 263 } else { 264 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 265 ierr = PetscMalloc1((n+1),&cia);CHKERRQ(ierr); 266 ierr = PetscMalloc1((nz+1),&cja);CHKERRQ(ierr); 267 jj = a->j; 268 for (i=0; i<nz; i++) { 269 collengths[jj[i]]++; 270 } 271 cia[0] = oshift; 272 for (i=0; i<n; i++) { 273 cia[i+1] = cia[i] + collengths[i]; 274 } 275 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 276 jj = a->j; 277 for (row=0; row<m; row++) { 278 mr = a->i[row+1] - a->i[row]; 279 for (i=0; i<mr; i++) { 280 col = *jj++; 281 282 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 283 } 284 } 285 ierr = PetscFree(collengths);CHKERRQ(ierr); 286 *ia = cia; *ja = cja; 287 } 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 293 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 294 { 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 if (!ia) PetscFunctionReturn(0); 299 300 ierr = PetscFree(*ia);CHKERRQ(ierr); 301 ierr = PetscFree(*ja);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 /* 306 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 307 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 308 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 309 */ 310 #undef __FUNCT__ 311 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 312 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 313 { 314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 315 PetscErrorCode ierr; 316 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 317 PetscInt nz = a->i[m],row,*jj,mr,col; 318 PetscInt *cspidx; 319 320 PetscFunctionBegin; 321 *nn = n; 322 if (!ia) PetscFunctionReturn(0); 323 324 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 325 ierr = PetscMalloc1((n+1),&cia);CHKERRQ(ierr); 326 ierr = PetscMalloc1((nz+1),&cja);CHKERRQ(ierr); 327 ierr = PetscMalloc1((nz+1),&cspidx);CHKERRQ(ierr); 328 jj = a->j; 329 for (i=0; i<nz; i++) { 330 collengths[jj[i]]++; 331 } 332 cia[0] = oshift; 333 for (i=0; i<n; i++) { 334 cia[i+1] = cia[i] + collengths[i]; 335 } 336 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 337 jj = a->j; 338 for (row=0; row<m; row++) { 339 mr = a->i[row+1] - a->i[row]; 340 for (i=0; i<mr; i++) { 341 col = *jj++; 342 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 343 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 344 } 345 } 346 ierr = PetscFree(collengths);CHKERRQ(ierr); 347 *ia = cia; *ja = cja; 348 *spidx = cspidx; 349 PetscFunctionReturn(0); 350 } 351 352 #undef __FUNCT__ 353 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 354 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 355 { 356 PetscErrorCode ierr; 357 358 PetscFunctionBegin; 359 ierr = MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);CHKERRQ(ierr); 360 ierr = PetscFree(*spidx);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 366 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 367 { 368 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 369 PetscInt *ai = a->i; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 374 PetscFunctionReturn(0); 375 } 376 377 #undef __FUNCT__ 378 #define __FUNCT__ "MatSeqAIJSetValuesLocalFast" 379 /* 380 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 381 382 - a single row of values is set with each call 383 - the local to global mapping is the identity 384 - no row or column indices are negative or (in error) larger than the number of rows or columns 385 - the values are always added to the matrix, not set 386 - the columns indices on each call are sorted 387 - no new locations are introduced in the nonzero structure of the matrix 388 389 */ 390 PetscErrorCode MatSeqAIJSetValuesLocalFastold(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 391 { 392 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 393 PetscInt low,high,t,row,nrow,i,col,l; 394 const PetscInt *rp,*ai = a->i,*ailen = a->ilen; 395 const PetscInt *aj = a->j; 396 MatScalar *ap,*aa = a->a; 397 PetscBool found; 398 399 PetscFunctionBegin; 400 if (m >1) printf("trouble1\n"); 401 row = im[0]; 402 rp = aj + ai[row]; 403 ap = aa + ai[row]; 404 nrow = ailen[row]; 405 printf("row %d nrow %d\n",row,nrow); 406 if (!nrow) printf("trouble4 %d\n",row); 407 low = 0; 408 for (l=0; l<n; l++) { /* loop over added columns */ 409 col = in[l]; 410 high = nrow; 411 found = PETSC_FALSE; 412 //printf("%d %d %d %d %d %d\n",l,low,high-1,rp[low],rp[high-1],col); 413 if (l > 0 && in[l] < in[l-1]) printf("trouble3 %d %d\n",in[l-1],in[l]); 414 while (high-low > 5) { 415 t = (low+high)/2; 416 if (rp[t] > col) high = t; 417 else low = t; 418 } 419 for (i=low; i<high; i++) { 420 if (rp[i] == col) { 421 ap[i] += v[l]; 422 low = i + 1; 423 found = PETSC_TRUE; 424 break; 425 } 426 } 427 if (!found) printf("trouble2 row %d %d %d %d %d %d %d %d\n",row,low,high-1,rp[low],rp[high-1],col,rp[0],rp[nrow-1]); 428 } 429 PetscFunctionReturn(0); 430 } 431 432 #include <petsc-private/isimpl.h> 433 #undef __FUNCT__ 434 #define __FUNCT__ "MatSeqAIJSetValuesLocalFast" 435 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 436 { 437 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 438 PetscInt *rp,low,high,t,ii,row,nrow,i,col,l,rmax,N; 439 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 440 PetscErrorCode ierr; 441 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 442 MatScalar *ap,value,*aa = a->a; 443 PetscBool ignorezeroentries = a->ignorezeroentries; 444 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 445 446 row = ridx[im[0]]; 447 rp = aj + ai[row]; ap = aa + ai[row]; 448 rmax = imax[row]; nrow = ailen[row]; 449 low = 0; 450 high = nrow; 451 for (l=0; l<n; l++) { /* loop over added columns */ 452 col = cidx[in[l]]; 453 value = v[l]; 454 455 if (col <= lastcol) low = 0; 456 else high = nrow; 457 lastcol = col; 458 while (high-low > 5) { 459 t = (low+high)/2; 460 if (rp[t] > col) high = t; 461 else low = t; 462 } 463 for (i=low; i<high; i++) { 464 if (rp[i] > col) break; 465 if (rp[i] == col) { 466 if (is == ADD_VALUES) ap[i] += value; 467 else ap[i] = value; 468 low = i + 1; 469 goto noinsert; 470 } 471 } 472 if (value == 0.0 && ignorezeroentries) goto noinsert; 473 if (nonew == 1) goto noinsert; 474 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 475 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 476 N = nrow++ - 1; a->nz++; high++; 477 /* shift up all the later entries in this row */ 478 for (ii=N; ii>=i; ii--) { 479 rp[ii+1] = rp[ii]; 480 ap[ii+1] = ap[ii]; 481 } 482 rp[i] = col; 483 ap[i] = value; 484 low = i + 1; 485 noinsert:; 486 } 487 return 0; 488 } 489 490 #undef __FUNCT__ 491 #define __FUNCT__ "MatSetValues_SeqAIJ" 492 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 493 { 494 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 495 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 496 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 497 PetscErrorCode ierr; 498 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 499 MatScalar *ap,value,*aa = a->a; 500 PetscBool ignorezeroentries = a->ignorezeroentries; 501 PetscBool roworiented = a->roworiented; 502 503 PetscFunctionBegin; 504 if (v) PetscValidScalarPointer(v,6); 505 for (k=0; k<m; k++) { /* loop over added rows */ 506 row = im[k]; 507 if (row < 0) continue; 508 #if defined(PETSC_USE_DEBUG) 509 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 510 #endif 511 rp = aj + ai[row]; ap = aa + ai[row]; 512 rmax = imax[row]; nrow = ailen[row]; 513 low = 0; 514 high = nrow; 515 for (l=0; l<n; l++) { /* loop over added columns */ 516 if (in[l] < 0) continue; 517 #if defined(PETSC_USE_DEBUG) 518 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 519 #endif 520 col = in[l]; 521 if (v) { 522 if (roworiented) { 523 value = v[l + k*n]; 524 } else { 525 value = v[k + l*m]; 526 } 527 } else { 528 value = 0.; 529 } 530 if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES)) continue; 531 532 if (col <= lastcol) low = 0; 533 else high = nrow; 534 lastcol = col; 535 while (high-low > 5) { 536 t = (low+high)/2; 537 if (rp[t] > col) high = t; 538 else low = t; 539 } 540 for (i=low; i<high; i++) { 541 if (rp[i] > col) break; 542 if (rp[i] == col) { 543 if (is == ADD_VALUES) ap[i] += value; 544 else ap[i] = value; 545 low = i + 1; 546 goto noinsert; 547 } 548 } 549 if (value == 0.0 && ignorezeroentries) goto noinsert; 550 if (nonew == 1) goto noinsert; 551 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 552 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 553 N = nrow++ - 1; a->nz++; high++; 554 /* shift up all the later entries in this row */ 555 for (ii=N; ii>=i; ii--) { 556 rp[ii+1] = rp[ii]; 557 ap[ii+1] = ap[ii]; 558 } 559 rp[i] = col; 560 ap[i] = value; 561 low = i + 1; 562 noinsert:; 563 } 564 ailen[row] = nrow; 565 } 566 A->same_nonzero = PETSC_FALSE; 567 PetscFunctionReturn(0); 568 } 569 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "MatGetValues_SeqAIJ" 573 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 574 { 575 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 576 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 577 PetscInt *ai = a->i,*ailen = a->ilen; 578 MatScalar *ap,*aa = a->a; 579 580 PetscFunctionBegin; 581 for (k=0; k<m; k++) { /* loop over rows */ 582 row = im[k]; 583 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 584 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 585 rp = aj + ai[row]; ap = aa + ai[row]; 586 nrow = ailen[row]; 587 for (l=0; l<n; l++) { /* loop over columns */ 588 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 589 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 590 col = in[l]; 591 high = nrow; low = 0; /* assume unsorted */ 592 while (high-low > 5) { 593 t = (low+high)/2; 594 if (rp[t] > col) high = t; 595 else low = t; 596 } 597 for (i=low; i<high; i++) { 598 if (rp[i] > col) break; 599 if (rp[i] == col) { 600 *v++ = ap[i]; 601 goto finished; 602 } 603 } 604 *v++ = 0.0; 605 finished:; 606 } 607 } 608 PetscFunctionReturn(0); 609 } 610 611 612 #undef __FUNCT__ 613 #define __FUNCT__ "MatView_SeqAIJ_Binary" 614 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 615 { 616 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 617 PetscErrorCode ierr; 618 PetscInt i,*col_lens; 619 int fd; 620 FILE *file; 621 622 PetscFunctionBegin; 623 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 624 ierr = PetscMalloc1((4+A->rmap->n),&col_lens);CHKERRQ(ierr); 625 626 col_lens[0] = MAT_FILE_CLASSID; 627 col_lens[1] = A->rmap->n; 628 col_lens[2] = A->cmap->n; 629 col_lens[3] = a->nz; 630 631 /* store lengths of each row and write (including header) to file */ 632 for (i=0; i<A->rmap->n; i++) { 633 col_lens[4+i] = a->i[i+1] - a->i[i]; 634 } 635 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 636 ierr = PetscFree(col_lens);CHKERRQ(ierr); 637 638 /* store column indices (zero start index) */ 639 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 640 641 /* store nonzero values */ 642 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 643 644 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 645 if (file) { 646 fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs); 647 } 648 PetscFunctionReturn(0); 649 } 650 651 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 655 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 656 { 657 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 658 PetscErrorCode ierr; 659 PetscInt i,j,m = A->rmap->n,shift=0; 660 const char *name; 661 PetscViewerFormat format; 662 663 PetscFunctionBegin; 664 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 665 if (format == PETSC_VIEWER_ASCII_MATLAB) { 666 PetscInt nofinalvalue = 0; 667 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift))) { 668 nofinalvalue = 1; 669 } 670 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 671 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 672 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 673 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 674 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 675 676 for (i=0; i<m; i++) { 677 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 678 #if defined(PETSC_USE_COMPLEX) 679 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 680 #else 681 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 682 #endif 683 } 684 } 685 if (nofinalvalue) { 686 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 687 } 688 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 689 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 690 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 691 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 692 PetscFunctionReturn(0); 693 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 694 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 695 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 696 for (i=0; i<m; i++) { 697 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 698 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 699 #if defined(PETSC_USE_COMPLEX) 700 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 701 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 702 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 703 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 704 } else if (PetscRealPart(a->a[j]) != 0.0) { 705 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 706 } 707 #else 708 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);CHKERRQ(ierr);} 709 #endif 710 } 711 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 712 } 713 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 714 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 715 PetscInt nzd=0,fshift=1,*sptr; 716 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 717 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 718 ierr = PetscMalloc1((m+1),&sptr);CHKERRQ(ierr); 719 for (i=0; i<m; i++) { 720 sptr[i] = nzd+1; 721 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 722 if (a->j[j] >= i) { 723 #if defined(PETSC_USE_COMPLEX) 724 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 725 #else 726 if (a->a[j] != 0.0) nzd++; 727 #endif 728 } 729 } 730 } 731 sptr[m] = nzd+1; 732 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 733 for (i=0; i<m+1; i+=6) { 734 if (i+4<m) { 735 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr); 736 } else if (i+3<m) { 737 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr); 738 } else if (i+2<m) { 739 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr); 740 } else if (i+1<m) { 741 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr); 742 } else if (i<m) { 743 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr); 744 } else { 745 ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr); 746 } 747 } 748 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 749 ierr = PetscFree(sptr);CHKERRQ(ierr); 750 for (i=0; i<m; i++) { 751 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 752 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 753 } 754 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 755 } 756 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 757 for (i=0; i<m; i++) { 758 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 759 if (a->j[j] >= i) { 760 #if defined(PETSC_USE_COMPLEX) 761 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 762 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 763 } 764 #else 765 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 766 #endif 767 } 768 } 769 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 770 } 771 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 772 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 773 PetscInt cnt = 0,jcnt; 774 PetscScalar value; 775 #if defined(PETSC_USE_COMPLEX) 776 PetscBool realonly = PETSC_TRUE; 777 778 for (i=0; i<a->i[m]; i++) { 779 if (PetscImaginaryPart(a->a[i]) != 0.0) { 780 realonly = PETSC_FALSE; 781 break; 782 } 783 } 784 #endif 785 786 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 787 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 788 for (i=0; i<m; i++) { 789 jcnt = 0; 790 for (j=0; j<A->cmap->n; j++) { 791 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 792 value = a->a[cnt++]; 793 jcnt++; 794 } else { 795 value = 0.0; 796 } 797 #if defined(PETSC_USE_COMPLEX) 798 if (realonly) { 799 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",PetscRealPart(value));CHKERRQ(ierr); 800 } else { 801 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 802 } 803 #else 804 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 805 #endif 806 } 807 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 808 } 809 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 810 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 811 PetscInt fshift=1; 812 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 813 #if defined(PETSC_USE_COMPLEX) 814 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 815 #else 816 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 817 #endif 818 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 819 for (i=0; i<m; i++) { 820 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 821 #if defined(PETSC_USE_COMPLEX) 822 if (PetscImaginaryPart(a->a[j]) > 0.0) { 823 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %g %g\n", i+fshift,a->j[j]+fshift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 824 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 825 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %g -%g\n", i+fshift,a->j[j]+fshift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 826 } else { 827 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %g\n", i+fshift,a->j[j]+fshift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 828 } 829 #else 830 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);CHKERRQ(ierr); 831 #endif 832 } 833 } 834 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 835 } else { 836 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 837 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 838 if (A->factortype) { 839 for (i=0; i<m; i++) { 840 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 841 /* L part */ 842 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 843 #if defined(PETSC_USE_COMPLEX) 844 if (PetscImaginaryPart(a->a[j]) > 0.0) { 845 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 846 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 847 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 848 } else { 849 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 850 } 851 #else 852 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);CHKERRQ(ierr); 853 #endif 854 } 855 /* diagonal */ 856 j = a->diag[i]; 857 #if defined(PETSC_USE_COMPLEX) 858 if (PetscImaginaryPart(a->a[j]) > 0.0) { 859 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),PetscImaginaryPart(1.0/a->a[j]));CHKERRQ(ierr); 860 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 861 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),-PetscImaginaryPart(1.0/a->a[j]));CHKERRQ(ierr); 862 } else { 863 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 864 } 865 #else 866 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)1.0/a->a[j]);CHKERRQ(ierr); 867 #endif 868 869 /* U part */ 870 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 871 #if defined(PETSC_USE_COMPLEX) 872 if (PetscImaginaryPart(a->a[j]) > 0.0) { 873 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 874 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 875 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 876 } else { 877 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 878 } 879 #else 880 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);CHKERRQ(ierr); 881 #endif 882 } 883 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 884 } 885 } else { 886 for (i=0; i<m; i++) { 887 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 888 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 889 #if defined(PETSC_USE_COMPLEX) 890 if (PetscImaginaryPart(a->a[j]) > 0.0) { 891 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 892 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 893 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 894 } else { 895 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 896 } 897 #else 898 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);CHKERRQ(ierr); 899 #endif 900 } 901 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 902 } 903 } 904 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 905 } 906 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 907 PetscFunctionReturn(0); 908 } 909 910 #include <petscdraw.h> 911 #undef __FUNCT__ 912 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 913 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 914 { 915 Mat A = (Mat) Aa; 916 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 917 PetscErrorCode ierr; 918 PetscInt i,j,m = A->rmap->n,color; 919 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 920 PetscViewer viewer; 921 PetscViewerFormat format; 922 923 PetscFunctionBegin; 924 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 925 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 926 927 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 928 /* loop over matrix elements drawing boxes */ 929 930 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 931 /* Blue for negative, Cyan for zero and Red for positive */ 932 color = PETSC_DRAW_BLUE; 933 for (i=0; i<m; i++) { 934 y_l = m - i - 1.0; y_r = y_l + 1.0; 935 for (j=a->i[i]; j<a->i[i+1]; j++) { 936 x_l = a->j[j]; x_r = x_l + 1.0; 937 if (PetscRealPart(a->a[j]) >= 0.) continue; 938 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 939 } 940 } 941 color = PETSC_DRAW_CYAN; 942 for (i=0; i<m; i++) { 943 y_l = m - i - 1.0; y_r = y_l + 1.0; 944 for (j=a->i[i]; j<a->i[i+1]; j++) { 945 x_l = a->j[j]; x_r = x_l + 1.0; 946 if (a->a[j] != 0.) continue; 947 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 948 } 949 } 950 color = PETSC_DRAW_RED; 951 for (i=0; i<m; i++) { 952 y_l = m - i - 1.0; y_r = y_l + 1.0; 953 for (j=a->i[i]; j<a->i[i+1]; j++) { 954 x_l = a->j[j]; x_r = x_l + 1.0; 955 if (PetscRealPart(a->a[j]) <= 0.) continue; 956 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 957 } 958 } 959 } else { 960 /* use contour shading to indicate magnitude of values */ 961 /* first determine max of all nonzero values */ 962 PetscInt nz = a->nz,count; 963 PetscDraw popup; 964 PetscReal scale; 965 966 for (i=0; i<nz; i++) { 967 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 968 } 969 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 970 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 971 if (popup) { 972 ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr); 973 } 974 count = 0; 975 for (i=0; i<m; i++) { 976 y_l = m - i - 1.0; y_r = y_l + 1.0; 977 for (j=a->i[i]; j<a->i[i+1]; j++) { 978 x_l = a->j[j]; x_r = x_l + 1.0; 979 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 980 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 981 count++; 982 } 983 } 984 } 985 PetscFunctionReturn(0); 986 } 987 988 #include <petscdraw.h> 989 #undef __FUNCT__ 990 #define __FUNCT__ "MatView_SeqAIJ_Draw" 991 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 992 { 993 PetscErrorCode ierr; 994 PetscDraw draw; 995 PetscReal xr,yr,xl,yl,h,w; 996 PetscBool isnull; 997 998 PetscFunctionBegin; 999 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1000 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1001 if (isnull) PetscFunctionReturn(0); 1002 1003 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1004 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1005 xr += w; yr += h; xl = -w; yl = -h; 1006 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1007 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 1008 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011 1012 #undef __FUNCT__ 1013 #define __FUNCT__ "MatView_SeqAIJ" 1014 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1015 { 1016 PetscErrorCode ierr; 1017 PetscBool iascii,isbinary,isdraw; 1018 1019 PetscFunctionBegin; 1020 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1021 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1022 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1023 if (iascii) { 1024 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 1025 } else if (isbinary) { 1026 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 1027 } else if (isdraw) { 1028 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 1029 } 1030 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 1036 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1037 { 1038 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1039 PetscErrorCode ierr; 1040 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 1041 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1042 MatScalar *aa = a->a,*ap; 1043 PetscReal ratio = 0.6; 1044 1045 PetscFunctionBegin; 1046 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1047 1048 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1049 for (i=1; i<m; i++) { 1050 /* move each row back by the amount of empty slots (fshift) before it*/ 1051 fshift += imax[i-1] - ailen[i-1]; 1052 rmax = PetscMax(rmax,ailen[i]); 1053 if (fshift) { 1054 ip = aj + ai[i]; 1055 ap = aa + ai[i]; 1056 N = ailen[i]; 1057 for (j=0; j<N; j++) { 1058 ip[j-fshift] = ip[j]; 1059 ap[j-fshift] = ap[j]; 1060 } 1061 } 1062 ai[i] = ai[i-1] + ailen[i-1]; 1063 } 1064 if (m) { 1065 fshift += imax[m-1] - ailen[m-1]; 1066 ai[m] = ai[m-1] + ailen[m-1]; 1067 } 1068 1069 /* reset ilen and imax for each row */ 1070 a->nonzerorowcnt = 0; 1071 for (i=0; i<m; i++) { 1072 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1073 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1074 } 1075 a->nz = ai[m]; 1076 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 1077 1078 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1079 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 1080 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 1081 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 1082 1083 A->info.mallocs += a->reallocs; 1084 a->reallocs = 0; 1085 A->info.nz_unneeded = (double)fshift; 1086 a->rmax = rmax; 1087 1088 ierr = MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 1089 1090 A->same_nonzero = PETSC_TRUE; 1091 1092 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 1093 1094 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1095 PetscFunctionReturn(0); 1096 } 1097 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "MatRealPart_SeqAIJ" 1100 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1101 { 1102 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1103 PetscInt i,nz = a->nz; 1104 MatScalar *aa = a->a; 1105 PetscErrorCode ierr; 1106 1107 PetscFunctionBegin; 1108 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1109 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 1115 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1116 { 1117 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1118 PetscInt i,nz = a->nz; 1119 MatScalar *aa = a->a; 1120 PetscErrorCode ierr; 1121 1122 PetscFunctionBegin; 1123 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1124 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 #if defined(PETSC_THREADCOMM_ACTIVE) 1129 PetscErrorCode MatZeroEntries_SeqAIJ_Kernel(PetscInt thread_id,Mat A) 1130 { 1131 PetscErrorCode ierr; 1132 PetscInt *trstarts=A->rmap->trstarts; 1133 PetscInt n,start,end; 1134 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1135 1136 start = trstarts[thread_id]; 1137 end = trstarts[thread_id+1]; 1138 n = a->i[end] - a->i[start]; 1139 ierr = PetscMemzero(a->a+a->i[start],n*sizeof(PetscScalar));CHKERRQ(ierr); 1140 return 0; 1141 } 1142 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1145 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1146 { 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatZeroEntries_SeqAIJ_Kernel,1,A);CHKERRQ(ierr); 1151 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1152 PetscFunctionReturn(0); 1153 } 1154 #else 1155 #undef __FUNCT__ 1156 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1157 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1158 { 1159 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1160 PetscErrorCode ierr; 1161 1162 PetscFunctionBegin; 1163 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1164 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 #endif 1168 1169 #undef __FUNCT__ 1170 #define __FUNCT__ "MatDestroy_SeqAIJ" 1171 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1172 { 1173 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1174 PetscErrorCode ierr; 1175 1176 PetscFunctionBegin; 1177 #if defined(PETSC_USE_LOG) 1178 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1179 #endif 1180 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1181 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1182 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1183 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1184 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1185 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1186 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1187 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1188 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1189 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1190 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1191 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 1192 ierr = MatDestroy(&a->XtoY);CHKERRQ(ierr); 1193 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1194 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1195 1196 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1197 ierr = PetscFree(A->data);CHKERRQ(ierr); 1198 1199 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1200 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1201 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1203 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1204 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1205 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1206 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1207 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1208 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1209 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "MatSetOption_SeqAIJ" 1215 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1216 { 1217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1218 PetscErrorCode ierr; 1219 1220 PetscFunctionBegin; 1221 switch (op) { 1222 case MAT_ROW_ORIENTED: 1223 a->roworiented = flg; 1224 break; 1225 case MAT_KEEP_NONZERO_PATTERN: 1226 a->keepnonzeropattern = flg; 1227 break; 1228 case MAT_NEW_NONZERO_LOCATIONS: 1229 a->nonew = (flg ? 0 : 1); 1230 break; 1231 case MAT_NEW_NONZERO_LOCATION_ERR: 1232 a->nonew = (flg ? -1 : 0); 1233 break; 1234 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1235 a->nonew = (flg ? -2 : 0); 1236 break; 1237 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1238 a->nounused = (flg ? -1 : 0); 1239 break; 1240 case MAT_IGNORE_ZERO_ENTRIES: 1241 a->ignorezeroentries = flg; 1242 break; 1243 case MAT_SPD: 1244 case MAT_SYMMETRIC: 1245 case MAT_STRUCTURALLY_SYMMETRIC: 1246 case MAT_HERMITIAN: 1247 case MAT_SYMMETRY_ETERNAL: 1248 /* These options are handled directly by MatSetOption() */ 1249 break; 1250 case MAT_NEW_DIAGONALS: 1251 case MAT_IGNORE_OFF_PROC_ENTRIES: 1252 case MAT_USE_HASH_TABLE: 1253 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1254 break; 1255 case MAT_USE_INODES: 1256 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1257 break; 1258 default: 1259 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1260 } 1261 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 1267 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1268 { 1269 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1270 PetscErrorCode ierr; 1271 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1272 PetscScalar *aa=a->a,*x,zero=0.0; 1273 1274 PetscFunctionBegin; 1275 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1276 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1277 1278 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1279 PetscInt *diag=a->diag; 1280 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1281 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1282 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 ierr = VecSet(v,zero);CHKERRQ(ierr); 1287 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1288 for (i=0; i<n; i++) { 1289 nz = ai[i+1] - ai[i]; 1290 if (!nz) x[i] = 0.0; 1291 for (j=ai[i]; j<ai[i+1]; j++) { 1292 if (aj[j] == i) { 1293 x[i] = aa[j]; 1294 break; 1295 } 1296 } 1297 } 1298 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 1305 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1306 { 1307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1308 PetscScalar *x,*y; 1309 PetscErrorCode ierr; 1310 PetscInt m = A->rmap->n; 1311 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1312 MatScalar *v; 1313 PetscScalar alpha; 1314 PetscInt n,i,j,*idx,*ii,*ridx=NULL; 1315 Mat_CompressedRow cprow = a->compressedrow; 1316 PetscBool usecprow = cprow.use; 1317 #endif 1318 1319 PetscFunctionBegin; 1320 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1321 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1322 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1323 1324 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1325 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1326 #else 1327 if (usecprow) { 1328 m = cprow.nrows; 1329 ii = cprow.i; 1330 ridx = cprow.rindex; 1331 } else { 1332 ii = a->i; 1333 } 1334 for (i=0; i<m; i++) { 1335 idx = a->j + ii[i]; 1336 v = a->a + ii[i]; 1337 n = ii[i+1] - ii[i]; 1338 if (usecprow) { 1339 alpha = x[ridx[i]]; 1340 } else { 1341 alpha = x[i]; 1342 } 1343 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1344 } 1345 #endif 1346 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1347 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1348 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1349 PetscFunctionReturn(0); 1350 } 1351 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 1354 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1355 { 1356 PetscErrorCode ierr; 1357 1358 PetscFunctionBegin; 1359 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1360 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1361 PetscFunctionReturn(0); 1362 } 1363 1364 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1365 #if defined(PETSC_THREADCOMM_ACTIVE) 1366 PetscErrorCode MatMult_SeqAIJ_Kernel(PetscInt thread_id,Mat A,Vec xx,Vec yy) 1367 { 1368 PetscErrorCode ierr; 1369 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1370 PetscScalar *y; 1371 const PetscScalar *x; 1372 const MatScalar *aa; 1373 PetscInt *trstarts=A->rmap->trstarts; 1374 PetscInt n,start,end,i; 1375 const PetscInt *aj,*ai; 1376 PetscScalar sum; 1377 1378 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1379 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1380 start = trstarts[thread_id]; 1381 end = trstarts[thread_id+1]; 1382 aj = a->j; 1383 aa = a->a; 1384 ai = a->i; 1385 for (i=start; i<end; i++) { 1386 n = ai[i+1] - ai[i]; 1387 aj = a->j + ai[i]; 1388 aa = a->a + ai[i]; 1389 sum = 0.0; 1390 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1391 y[i] = sum; 1392 } 1393 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1394 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1395 return 0; 1396 } 1397 1398 #undef __FUNCT__ 1399 #define __FUNCT__ "MatMult_SeqAIJ" 1400 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1401 { 1402 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1403 PetscScalar *y; 1404 const PetscScalar *x; 1405 const MatScalar *aa; 1406 PetscErrorCode ierr; 1407 PetscInt m=A->rmap->n; 1408 const PetscInt *aj,*ii,*ridx=NULL; 1409 PetscInt n,i; 1410 PetscScalar sum; 1411 PetscBool usecprow=a->compressedrow.use; 1412 1413 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1414 #pragma disjoint(*x,*y,*aa) 1415 #endif 1416 1417 PetscFunctionBegin; 1418 aj = a->j; 1419 aa = a->a; 1420 ii = a->i; 1421 if (usecprow) { /* use compressed row format */ 1422 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1423 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1424 m = a->compressedrow.nrows; 1425 ii = a->compressedrow.i; 1426 ridx = a->compressedrow.rindex; 1427 for (i=0; i<m; i++) { 1428 n = ii[i+1] - ii[i]; 1429 aj = a->j + ii[i]; 1430 aa = a->a + ii[i]; 1431 sum = 0.0; 1432 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1433 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1434 y[*ridx++] = sum; 1435 } 1436 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1437 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1438 } else { /* do not use compressed row format */ 1439 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1440 fortranmultaij_(&m,x,ii,aj,aa,y); 1441 #else 1442 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1443 #endif 1444 } 1445 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1446 PetscFunctionReturn(0); 1447 } 1448 #else 1449 #undef __FUNCT__ 1450 #define __FUNCT__ "MatMult_SeqAIJ" 1451 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1452 { 1453 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1454 PetscScalar *y; 1455 const PetscScalar *x; 1456 const MatScalar *aa; 1457 PetscErrorCode ierr; 1458 PetscInt m=A->rmap->n; 1459 const PetscInt *aj,*ii,*ridx=NULL; 1460 PetscInt n,i; 1461 PetscScalar sum; 1462 PetscBool usecprow=a->compressedrow.use; 1463 1464 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1465 #pragma disjoint(*x,*y,*aa) 1466 #endif 1467 1468 PetscFunctionBegin; 1469 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1470 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1471 aj = a->j; 1472 aa = a->a; 1473 ii = a->i; 1474 if (usecprow) { /* use compressed row format */ 1475 m = a->compressedrow.nrows; 1476 ii = a->compressedrow.i; 1477 ridx = a->compressedrow.rindex; 1478 for (i=0; i<m; i++) { 1479 n = ii[i+1] - ii[i]; 1480 aj = a->j + ii[i]; 1481 aa = a->a + ii[i]; 1482 sum = 0.0; 1483 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1484 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1485 y[*ridx++] = sum; 1486 } 1487 } else { /* do not use compressed row format */ 1488 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1489 fortranmultaij_(&m,x,ii,aj,aa,y); 1490 #else 1491 #if defined(PETSC_THREADCOMM_ACTIVE) 1492 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1493 #else 1494 for (i=0; i<m; i++) { 1495 n = ii[i+1] - ii[i]; 1496 aj = a->j + ii[i]; 1497 aa = a->a + ii[i]; 1498 sum = 0.0; 1499 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1500 y[i] = sum; 1501 } 1502 #endif 1503 #endif 1504 } 1505 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1506 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1507 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 #endif 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "MatMultMax_SeqAIJ" 1514 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1515 { 1516 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1517 PetscScalar *y; 1518 const PetscScalar *x; 1519 const MatScalar *aa; 1520 PetscErrorCode ierr; 1521 PetscInt m=A->rmap->n; 1522 const PetscInt *aj,*ii,*ridx=NULL; 1523 PetscInt n,i,nonzerorow=0; 1524 PetscScalar sum; 1525 PetscBool usecprow=a->compressedrow.use; 1526 1527 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1528 #pragma disjoint(*x,*y,*aa) 1529 #endif 1530 1531 PetscFunctionBegin; 1532 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1533 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1534 aj = a->j; 1535 aa = a->a; 1536 ii = a->i; 1537 if (usecprow) { /* use compressed row format */ 1538 m = a->compressedrow.nrows; 1539 ii = a->compressedrow.i; 1540 ridx = a->compressedrow.rindex; 1541 for (i=0; i<m; i++) { 1542 n = ii[i+1] - ii[i]; 1543 aj = a->j + ii[i]; 1544 aa = a->a + ii[i]; 1545 sum = 0.0; 1546 nonzerorow += (n>0); 1547 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1548 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1549 y[*ridx++] = sum; 1550 } 1551 } else { /* do not use compressed row format */ 1552 for (i=0; i<m; i++) { 1553 n = ii[i+1] - ii[i]; 1554 aj = a->j + ii[i]; 1555 aa = a->a + ii[i]; 1556 sum = 0.0; 1557 nonzerorow += (n>0); 1558 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1559 y[i] = sum; 1560 } 1561 } 1562 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1563 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1564 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #undef __FUNCT__ 1569 #define __FUNCT__ "MatMultAddMax_SeqAIJ" 1570 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1571 { 1572 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1573 PetscScalar *y,*z; 1574 const PetscScalar *x; 1575 const MatScalar *aa; 1576 PetscErrorCode ierr; 1577 PetscInt m = A->rmap->n,*aj,*ii; 1578 PetscInt n,i,*ridx=NULL; 1579 PetscScalar sum; 1580 PetscBool usecprow=a->compressedrow.use; 1581 1582 PetscFunctionBegin; 1583 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1584 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1585 if (zz != yy) { 1586 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1587 } else { 1588 z = y; 1589 } 1590 1591 aj = a->j; 1592 aa = a->a; 1593 ii = a->i; 1594 if (usecprow) { /* use compressed row format */ 1595 if (zz != yy) { 1596 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1597 } 1598 m = a->compressedrow.nrows; 1599 ii = a->compressedrow.i; 1600 ridx = a->compressedrow.rindex; 1601 for (i=0; i<m; i++) { 1602 n = ii[i+1] - ii[i]; 1603 aj = a->j + ii[i]; 1604 aa = a->a + ii[i]; 1605 sum = y[*ridx]; 1606 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1607 z[*ridx++] = sum; 1608 } 1609 } else { /* do not use compressed row format */ 1610 for (i=0; i<m; i++) { 1611 n = ii[i+1] - ii[i]; 1612 aj = a->j + ii[i]; 1613 aa = a->a + ii[i]; 1614 sum = y[i]; 1615 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1616 z[i] = sum; 1617 } 1618 } 1619 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1620 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1621 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1622 if (zz != yy) { 1623 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1624 } 1625 PetscFunctionReturn(0); 1626 } 1627 1628 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1629 #undef __FUNCT__ 1630 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1631 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1632 { 1633 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1634 PetscScalar *y,*z; 1635 const PetscScalar *x; 1636 const MatScalar *aa; 1637 PetscErrorCode ierr; 1638 PetscInt m = A->rmap->n,*aj,*ii; 1639 PetscInt n,i,*ridx=NULL; 1640 PetscScalar sum; 1641 PetscBool usecprow=a->compressedrow.use; 1642 1643 PetscFunctionBegin; 1644 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1645 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1646 if (zz != yy) { 1647 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1648 } else { 1649 z = y; 1650 } 1651 1652 aj = a->j; 1653 aa = a->a; 1654 ii = a->i; 1655 if (usecprow) { /* use compressed row format */ 1656 if (zz != yy) { 1657 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1658 } 1659 m = a->compressedrow.nrows; 1660 ii = a->compressedrow.i; 1661 ridx = a->compressedrow.rindex; 1662 for (i=0; i<m; i++) { 1663 n = ii[i+1] - ii[i]; 1664 aj = a->j + ii[i]; 1665 aa = a->a + ii[i]; 1666 sum = y[*ridx]; 1667 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1668 z[*ridx++] = sum; 1669 } 1670 } else { /* do not use compressed row format */ 1671 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1672 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1673 #else 1674 for (i=0; i<m; i++) { 1675 n = ii[i+1] - ii[i]; 1676 aj = a->j + ii[i]; 1677 aa = a->a + ii[i]; 1678 sum = y[i]; 1679 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1680 z[i] = sum; 1681 } 1682 #endif 1683 } 1684 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1685 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1686 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1687 if (zz != yy) { 1688 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1689 } 1690 #if defined(PETSC_HAVE_CUSP) 1691 /* 1692 ierr = VecView(xx,0);CHKERRQ(ierr); 1693 ierr = VecView(zz,0);CHKERRQ(ierr); 1694 ierr = MatView(A,0);CHKERRQ(ierr); 1695 */ 1696 #endif 1697 PetscFunctionReturn(0); 1698 } 1699 1700 /* 1701 Adds diagonal pointers to sparse matrix structure. 1702 */ 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1705 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1706 { 1707 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1708 PetscErrorCode ierr; 1709 PetscInt i,j,m = A->rmap->n; 1710 1711 PetscFunctionBegin; 1712 if (!a->diag) { 1713 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1714 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1715 } 1716 for (i=0; i<A->rmap->n; i++) { 1717 a->diag[i] = a->i[i+1]; 1718 for (j=a->i[i]; j<a->i[i+1]; j++) { 1719 if (a->j[j] == i) { 1720 a->diag[i] = j; 1721 break; 1722 } 1723 } 1724 } 1725 PetscFunctionReturn(0); 1726 } 1727 1728 /* 1729 Checks for missing diagonals 1730 */ 1731 #undef __FUNCT__ 1732 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1733 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1734 { 1735 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1736 PetscInt *diag,*jj = a->j,i; 1737 1738 PetscFunctionBegin; 1739 *missing = PETSC_FALSE; 1740 if (A->rmap->n > 0 && !jj) { 1741 *missing = PETSC_TRUE; 1742 if (d) *d = 0; 1743 PetscInfo(A,"Matrix has no entries therefore is missing diagonal"); 1744 } else { 1745 diag = a->diag; 1746 for (i=0; i<A->rmap->n; i++) { 1747 if (jj[diag[i]] != i) { 1748 *missing = PETSC_TRUE; 1749 if (d) *d = i; 1750 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1751 break; 1752 } 1753 } 1754 } 1755 PetscFunctionReturn(0); 1756 } 1757 1758 #undef __FUNCT__ 1759 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1760 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1761 { 1762 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1763 PetscErrorCode ierr; 1764 PetscInt i,*diag,m = A->rmap->n; 1765 MatScalar *v = a->a; 1766 PetscScalar *idiag,*mdiag; 1767 1768 PetscFunctionBegin; 1769 if (a->idiagvalid) PetscFunctionReturn(0); 1770 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1771 diag = a->diag; 1772 if (!a->idiag) { 1773 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1774 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1775 v = a->a; 1776 } 1777 mdiag = a->mdiag; 1778 idiag = a->idiag; 1779 1780 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1781 for (i=0; i<m; i++) { 1782 mdiag[i] = v[diag[i]]; 1783 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1784 idiag[i] = 1.0/v[diag[i]]; 1785 } 1786 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1787 } else { 1788 for (i=0; i<m; i++) { 1789 mdiag[i] = v[diag[i]]; 1790 idiag[i] = omega/(fshift + v[diag[i]]); 1791 } 1792 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1793 } 1794 a->idiagvalid = PETSC_TRUE; 1795 PetscFunctionReturn(0); 1796 } 1797 1798 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1799 #undef __FUNCT__ 1800 #define __FUNCT__ "MatSOR_SeqAIJ" 1801 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1802 { 1803 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1804 PetscScalar *x,d,sum,*t,scale; 1805 const MatScalar *v = a->a,*idiag=0,*mdiag; 1806 const PetscScalar *b, *bs,*xb, *ts; 1807 PetscErrorCode ierr; 1808 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1809 const PetscInt *idx,*diag; 1810 1811 PetscFunctionBegin; 1812 its = its*lits; 1813 1814 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1815 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1816 a->fshift = fshift; 1817 a->omega = omega; 1818 1819 diag = a->diag; 1820 t = a->ssor_work; 1821 idiag = a->idiag; 1822 mdiag = a->mdiag; 1823 1824 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1825 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1826 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1827 if (flag == SOR_APPLY_UPPER) { 1828 /* apply (U + D/omega) to the vector */ 1829 bs = b; 1830 for (i=0; i<m; i++) { 1831 d = fshift + mdiag[i]; 1832 n = a->i[i+1] - diag[i] - 1; 1833 idx = a->j + diag[i] + 1; 1834 v = a->a + diag[i] + 1; 1835 sum = b[i]*d/omega; 1836 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1837 x[i] = sum; 1838 } 1839 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1840 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1841 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1846 else if (flag & SOR_EISENSTAT) { 1847 /* Let A = L + U + D; where L is lower trianglar, 1848 U is upper triangular, E = D/omega; This routine applies 1849 1850 (L + E)^{-1} A (U + E)^{-1} 1851 1852 to a vector efficiently using Eisenstat's trick. 1853 */ 1854 scale = (2.0/omega) - 1.0; 1855 1856 /* x = (E + U)^{-1} b */ 1857 for (i=m-1; i>=0; i--) { 1858 n = a->i[i+1] - diag[i] - 1; 1859 idx = a->j + diag[i] + 1; 1860 v = a->a + diag[i] + 1; 1861 sum = b[i]; 1862 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1863 x[i] = sum*idiag[i]; 1864 } 1865 1866 /* t = b - (2*E - D)x */ 1867 v = a->a; 1868 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1869 1870 /* t = (E + L)^{-1}t */ 1871 ts = t; 1872 diag = a->diag; 1873 for (i=0; i<m; i++) { 1874 n = diag[i] - a->i[i]; 1875 idx = a->j + a->i[i]; 1876 v = a->a + a->i[i]; 1877 sum = t[i]; 1878 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1879 t[i] = sum*idiag[i]; 1880 /* x = x + t */ 1881 x[i] += t[i]; 1882 } 1883 1884 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1885 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1886 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1887 PetscFunctionReturn(0); 1888 } 1889 if (flag & SOR_ZERO_INITIAL_GUESS) { 1890 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1891 for (i=0; i<m; i++) { 1892 n = diag[i] - a->i[i]; 1893 idx = a->j + a->i[i]; 1894 v = a->a + a->i[i]; 1895 sum = b[i]; 1896 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1897 t[i] = sum; 1898 x[i] = sum*idiag[i]; 1899 } 1900 xb = t; 1901 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1902 } else xb = b; 1903 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1904 for (i=m-1; i>=0; i--) { 1905 n = a->i[i+1] - diag[i] - 1; 1906 idx = a->j + diag[i] + 1; 1907 v = a->a + diag[i] + 1; 1908 sum = xb[i]; 1909 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1910 if (xb == b) { 1911 x[i] = sum*idiag[i]; 1912 } else { 1913 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1914 } 1915 } 1916 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1917 } 1918 its--; 1919 } 1920 while (its--) { 1921 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1922 for (i=0; i<m; i++) { 1923 /* lower */ 1924 n = diag[i] - a->i[i]; 1925 idx = a->j + a->i[i]; 1926 v = a->a + a->i[i]; 1927 sum = b[i]; 1928 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1929 t[i] = sum; /* save application of the lower-triangular part */ 1930 /* upper */ 1931 n = a->i[i+1] - diag[i] - 1; 1932 idx = a->j + diag[i] + 1; 1933 v = a->a + diag[i] + 1; 1934 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1935 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1936 } 1937 xb = t; 1938 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1939 } else xb = b; 1940 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1941 for (i=m-1; i>=0; i--) { 1942 sum = xb[i]; 1943 if (xb == b) { 1944 /* whole matrix (no checkpointing available) */ 1945 n = a->i[i+1] - a->i[i]; 1946 idx = a->j + a->i[i]; 1947 v = a->a + a->i[i]; 1948 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1949 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1950 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1951 n = a->i[i+1] - diag[i] - 1; 1952 idx = a->j + diag[i] + 1; 1953 v = a->a + diag[i] + 1; 1954 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1955 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1956 } 1957 } 1958 if (xb == b) { 1959 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1960 } else { 1961 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1962 } 1963 } 1964 } 1965 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1966 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1967 PetscFunctionReturn(0); 1968 } 1969 1970 1971 #undef __FUNCT__ 1972 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1973 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1974 { 1975 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1976 1977 PetscFunctionBegin; 1978 info->block_size = 1.0; 1979 info->nz_allocated = (double)a->maxnz; 1980 info->nz_used = (double)a->nz; 1981 info->nz_unneeded = (double)(a->maxnz - a->nz); 1982 info->assemblies = (double)A->num_ass; 1983 info->mallocs = (double)A->info.mallocs; 1984 info->memory = ((PetscObject)A)->mem; 1985 if (A->factortype) { 1986 info->fill_ratio_given = A->info.fill_ratio_given; 1987 info->fill_ratio_needed = A->info.fill_ratio_needed; 1988 info->factor_mallocs = A->info.factor_mallocs; 1989 } else { 1990 info->fill_ratio_given = 0; 1991 info->fill_ratio_needed = 0; 1992 info->factor_mallocs = 0; 1993 } 1994 PetscFunctionReturn(0); 1995 } 1996 1997 #undef __FUNCT__ 1998 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1999 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2000 { 2001 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2002 PetscInt i,m = A->rmap->n - 1,d = 0; 2003 PetscErrorCode ierr; 2004 const PetscScalar *xx; 2005 PetscScalar *bb; 2006 PetscBool missing; 2007 2008 PetscFunctionBegin; 2009 if (x && b) { 2010 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2011 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2012 for (i=0; i<N; i++) { 2013 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2014 bb[rows[i]] = diag*xx[rows[i]]; 2015 } 2016 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2017 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2018 } 2019 2020 if (a->keepnonzeropattern) { 2021 for (i=0; i<N; i++) { 2022 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2023 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2024 } 2025 if (diag != 0.0) { 2026 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2027 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 2028 for (i=0; i<N; i++) { 2029 a->a[a->diag[rows[i]]] = diag; 2030 } 2031 } 2032 A->same_nonzero = PETSC_TRUE; 2033 } else { 2034 if (diag != 0.0) { 2035 for (i=0; i<N; i++) { 2036 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2037 if (a->ilen[rows[i]] > 0) { 2038 a->ilen[rows[i]] = 1; 2039 a->a[a->i[rows[i]]] = diag; 2040 a->j[a->i[rows[i]]] = rows[i]; 2041 } else { /* in case row was completely empty */ 2042 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 2043 } 2044 } 2045 } else { 2046 for (i=0; i<N; i++) { 2047 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2048 a->ilen[rows[i]] = 0; 2049 } 2050 } 2051 A->same_nonzero = PETSC_FALSE; 2052 } 2053 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 2057 #undef __FUNCT__ 2058 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 2059 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2060 { 2061 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2062 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2063 PetscErrorCode ierr; 2064 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2065 const PetscScalar *xx; 2066 PetscScalar *bb; 2067 2068 PetscFunctionBegin; 2069 if (x && b) { 2070 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2071 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2072 vecs = PETSC_TRUE; 2073 } 2074 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2075 for (i=0; i<N; i++) { 2076 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2077 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2078 2079 zeroed[rows[i]] = PETSC_TRUE; 2080 } 2081 for (i=0; i<A->rmap->n; i++) { 2082 if (!zeroed[i]) { 2083 for (j=a->i[i]; j<a->i[i+1]; j++) { 2084 if (zeroed[a->j[j]]) { 2085 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2086 a->a[j] = 0.0; 2087 } 2088 } 2089 } else if (vecs) bb[i] = diag*xx[i]; 2090 } 2091 if (x && b) { 2092 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2093 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2094 } 2095 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2096 if (diag != 0.0) { 2097 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2098 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 2099 for (i=0; i<N; i++) { 2100 a->a[a->diag[rows[i]]] = diag; 2101 } 2102 } 2103 A->same_nonzero = PETSC_TRUE; 2104 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "MatGetRow_SeqAIJ" 2110 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2111 { 2112 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2113 PetscInt *itmp; 2114 2115 PetscFunctionBegin; 2116 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2117 2118 *nz = a->i[row+1] - a->i[row]; 2119 if (v) *v = a->a + a->i[row]; 2120 if (idx) { 2121 itmp = a->j + a->i[row]; 2122 if (*nz) *idx = itmp; 2123 else *idx = 0; 2124 } 2125 PetscFunctionReturn(0); 2126 } 2127 2128 /* remove this function? */ 2129 #undef __FUNCT__ 2130 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 2131 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2132 { 2133 PetscFunctionBegin; 2134 PetscFunctionReturn(0); 2135 } 2136 2137 #undef __FUNCT__ 2138 #define __FUNCT__ "MatNorm_SeqAIJ" 2139 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2140 { 2141 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2142 MatScalar *v = a->a; 2143 PetscReal sum = 0.0; 2144 PetscErrorCode ierr; 2145 PetscInt i,j; 2146 2147 PetscFunctionBegin; 2148 if (type == NORM_FROBENIUS) { 2149 for (i=0; i<a->nz; i++) { 2150 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2151 } 2152 *nrm = PetscSqrtReal(sum); 2153 } else if (type == NORM_1) { 2154 PetscReal *tmp; 2155 PetscInt *jj = a->j; 2156 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2157 *nrm = 0.0; 2158 for (j=0; j<a->nz; j++) { 2159 tmp[*jj++] += PetscAbsScalar(*v); v++; 2160 } 2161 for (j=0; j<A->cmap->n; j++) { 2162 if (tmp[j] > *nrm) *nrm = tmp[j]; 2163 } 2164 ierr = PetscFree(tmp);CHKERRQ(ierr); 2165 } else if (type == NORM_INFINITY) { 2166 *nrm = 0.0; 2167 for (j=0; j<A->rmap->n; j++) { 2168 v = a->a + a->i[j]; 2169 sum = 0.0; 2170 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2171 sum += PetscAbsScalar(*v); v++; 2172 } 2173 if (sum > *nrm) *nrm = sum; 2174 } 2175 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2176 PetscFunctionReturn(0); 2177 } 2178 2179 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2180 #undef __FUNCT__ 2181 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 2182 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2183 { 2184 PetscErrorCode ierr; 2185 PetscInt i,j,anzj; 2186 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2187 PetscInt an=A->cmap->N,am=A->rmap->N; 2188 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2189 2190 PetscFunctionBegin; 2191 /* Allocate space for symbolic transpose info and work array */ 2192 ierr = PetscCalloc1((an+1),&ati);CHKERRQ(ierr); 2193 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2194 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2195 2196 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2197 /* Note: offset by 1 for fast conversion into csr format. */ 2198 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2199 /* Form ati for csr format of A^T. */ 2200 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2201 2202 /* Copy ati into atfill so we have locations of the next free space in atj */ 2203 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2204 2205 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2206 for (i=0;i<am;i++) { 2207 anzj = ai[i+1] - ai[i]; 2208 for (j=0;j<anzj;j++) { 2209 atj[atfill[*aj]] = i; 2210 atfill[*aj++] += 1; 2211 } 2212 } 2213 2214 /* Clean up temporary space and complete requests. */ 2215 ierr = PetscFree(atfill);CHKERRQ(ierr); 2216 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2217 2218 (*B)->rmap->bs = A->cmap->bs; 2219 (*B)->cmap->bs = A->rmap->bs; 2220 2221 b = (Mat_SeqAIJ*)((*B)->data); 2222 b->free_a = PETSC_FALSE; 2223 b->free_ij = PETSC_TRUE; 2224 b->nonew = 0; 2225 PetscFunctionReturn(0); 2226 } 2227 2228 #undef __FUNCT__ 2229 #define __FUNCT__ "MatTranspose_SeqAIJ" 2230 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2231 { 2232 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2233 Mat C; 2234 PetscErrorCode ierr; 2235 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2236 MatScalar *array = a->a; 2237 2238 PetscFunctionBegin; 2239 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 2240 2241 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2242 ierr = PetscCalloc1((1+A->cmap->n),&col);CHKERRQ(ierr); 2243 2244 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2245 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2246 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2247 ierr = MatSetBlockSizes(C,A->cmap->bs,A->rmap->bs);CHKERRQ(ierr); 2248 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2249 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2250 ierr = PetscFree(col);CHKERRQ(ierr); 2251 } else { 2252 C = *B; 2253 } 2254 2255 for (i=0; i<m; i++) { 2256 len = ai[i+1]-ai[i]; 2257 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2258 array += len; 2259 aj += len; 2260 } 2261 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2262 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2263 2264 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2265 *B = C; 2266 } else { 2267 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 2268 } 2269 PetscFunctionReturn(0); 2270 } 2271 2272 #undef __FUNCT__ 2273 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2274 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2275 { 2276 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2277 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2278 MatScalar *va,*vb; 2279 PetscErrorCode ierr; 2280 PetscInt ma,na,mb,nb, i; 2281 2282 PetscFunctionBegin; 2283 bij = (Mat_SeqAIJ*) B->data; 2284 2285 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2286 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2287 if (ma!=nb || na!=mb) { 2288 *f = PETSC_FALSE; 2289 PetscFunctionReturn(0); 2290 } 2291 aii = aij->i; bii = bij->i; 2292 adx = aij->j; bdx = bij->j; 2293 va = aij->a; vb = bij->a; 2294 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2295 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2296 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2297 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2298 2299 *f = PETSC_TRUE; 2300 for (i=0; i<ma; i++) { 2301 while (aptr[i]<aii[i+1]) { 2302 PetscInt idc,idr; 2303 PetscScalar vc,vr; 2304 /* column/row index/value */ 2305 idc = adx[aptr[i]]; 2306 idr = bdx[bptr[idc]]; 2307 vc = va[aptr[i]]; 2308 vr = vb[bptr[idc]]; 2309 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2310 *f = PETSC_FALSE; 2311 goto done; 2312 } else { 2313 aptr[i]++; 2314 if (B || i!=idc) bptr[idc]++; 2315 } 2316 } 2317 } 2318 done: 2319 ierr = PetscFree(aptr);CHKERRQ(ierr); 2320 ierr = PetscFree(bptr);CHKERRQ(ierr); 2321 PetscFunctionReturn(0); 2322 } 2323 2324 #undef __FUNCT__ 2325 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2326 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2327 { 2328 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2329 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2330 MatScalar *va,*vb; 2331 PetscErrorCode ierr; 2332 PetscInt ma,na,mb,nb, i; 2333 2334 PetscFunctionBegin; 2335 bij = (Mat_SeqAIJ*) B->data; 2336 2337 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2338 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2339 if (ma!=nb || na!=mb) { 2340 *f = PETSC_FALSE; 2341 PetscFunctionReturn(0); 2342 } 2343 aii = aij->i; bii = bij->i; 2344 adx = aij->j; bdx = bij->j; 2345 va = aij->a; vb = bij->a; 2346 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2347 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2348 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2349 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2350 2351 *f = PETSC_TRUE; 2352 for (i=0; i<ma; i++) { 2353 while (aptr[i]<aii[i+1]) { 2354 PetscInt idc,idr; 2355 PetscScalar vc,vr; 2356 /* column/row index/value */ 2357 idc = adx[aptr[i]]; 2358 idr = bdx[bptr[idc]]; 2359 vc = va[aptr[i]]; 2360 vr = vb[bptr[idc]]; 2361 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2362 *f = PETSC_FALSE; 2363 goto done; 2364 } else { 2365 aptr[i]++; 2366 if (B || i!=idc) bptr[idc]++; 2367 } 2368 } 2369 } 2370 done: 2371 ierr = PetscFree(aptr);CHKERRQ(ierr); 2372 ierr = PetscFree(bptr);CHKERRQ(ierr); 2373 PetscFunctionReturn(0); 2374 } 2375 2376 #undef __FUNCT__ 2377 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2378 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2379 { 2380 PetscErrorCode ierr; 2381 2382 PetscFunctionBegin; 2383 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2384 PetscFunctionReturn(0); 2385 } 2386 2387 #undef __FUNCT__ 2388 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2389 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2390 { 2391 PetscErrorCode ierr; 2392 2393 PetscFunctionBegin; 2394 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2395 PetscFunctionReturn(0); 2396 } 2397 2398 #undef __FUNCT__ 2399 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2400 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2401 { 2402 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2403 PetscScalar *l,*r,x; 2404 MatScalar *v; 2405 PetscErrorCode ierr; 2406 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2407 2408 PetscFunctionBegin; 2409 if (ll) { 2410 /* The local size is used so that VecMPI can be passed to this routine 2411 by MatDiagonalScale_MPIAIJ */ 2412 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2413 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2414 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2415 v = a->a; 2416 for (i=0; i<m; i++) { 2417 x = l[i]; 2418 M = a->i[i+1] - a->i[i]; 2419 for (j=0; j<M; j++) (*v++) *= x; 2420 } 2421 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2422 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2423 } 2424 if (rr) { 2425 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2426 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2427 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2428 v = a->a; jj = a->j; 2429 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2430 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2431 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2432 } 2433 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2434 PetscFunctionReturn(0); 2435 } 2436 2437 #undef __FUNCT__ 2438 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2439 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2440 { 2441 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2442 PetscErrorCode ierr; 2443 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2444 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2445 const PetscInt *irow,*icol; 2446 PetscInt nrows,ncols; 2447 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2448 MatScalar *a_new,*mat_a; 2449 Mat C; 2450 PetscBool stride,sorted; 2451 2452 PetscFunctionBegin; 2453 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2454 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2455 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2456 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2457 2458 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2459 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2460 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2461 2462 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2463 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2464 if (stride && step == 1) { 2465 /* special case of contiguous rows */ 2466 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2467 /* loop over new rows determining lens and starting points */ 2468 for (i=0; i<nrows; i++) { 2469 kstart = ai[irow[i]]; 2470 kend = kstart + ailen[irow[i]]; 2471 for (k=kstart; k<kend; k++) { 2472 if (aj[k] >= first) { 2473 starts[i] = k; 2474 break; 2475 } 2476 } 2477 sum = 0; 2478 while (k < kend) { 2479 if (aj[k++] >= first+ncols) break; 2480 sum++; 2481 } 2482 lens[i] = sum; 2483 } 2484 /* create submatrix */ 2485 if (scall == MAT_REUSE_MATRIX) { 2486 PetscInt n_cols,n_rows; 2487 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2488 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2489 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2490 C = *B; 2491 } else { 2492 PetscInt rbs,cbs; 2493 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2494 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2495 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2496 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2497 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2498 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2499 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2500 } 2501 c = (Mat_SeqAIJ*)C->data; 2502 2503 /* loop over rows inserting into submatrix */ 2504 a_new = c->a; 2505 j_new = c->j; 2506 i_new = c->i; 2507 2508 for (i=0; i<nrows; i++) { 2509 ii = starts[i]; 2510 lensi = lens[i]; 2511 for (k=0; k<lensi; k++) { 2512 *j_new++ = aj[ii+k] - first; 2513 } 2514 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2515 a_new += lensi; 2516 i_new[i+1] = i_new[i] + lensi; 2517 c->ilen[i] = lensi; 2518 } 2519 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2520 } else { 2521 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2522 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2523 ierr = PetscMalloc1((1+nrows),&lens);CHKERRQ(ierr); 2524 for (i=0; i<ncols; i++) { 2525 #if defined(PETSC_USE_DEBUG) 2526 if (icol[i] >= oldcols) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%D] %D <= A->cmap->n %D",i,icol[i],oldcols); 2527 #endif 2528 smap[icol[i]] = i+1; 2529 } 2530 2531 /* determine lens of each row */ 2532 for (i=0; i<nrows; i++) { 2533 kstart = ai[irow[i]]; 2534 kend = kstart + a->ilen[irow[i]]; 2535 lens[i] = 0; 2536 for (k=kstart; k<kend; k++) { 2537 if (smap[aj[k]]) { 2538 lens[i]++; 2539 } 2540 } 2541 } 2542 /* Create and fill new matrix */ 2543 if (scall == MAT_REUSE_MATRIX) { 2544 PetscBool equal; 2545 2546 c = (Mat_SeqAIJ*)((*B)->data); 2547 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2548 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2549 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2550 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2551 C = *B; 2552 } else { 2553 PetscInt rbs,cbs; 2554 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2555 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2556 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2557 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2558 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2559 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2560 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2561 } 2562 c = (Mat_SeqAIJ*)(C->data); 2563 for (i=0; i<nrows; i++) { 2564 row = irow[i]; 2565 kstart = ai[row]; 2566 kend = kstart + a->ilen[row]; 2567 mat_i = c->i[i]; 2568 mat_j = c->j + mat_i; 2569 mat_a = c->a + mat_i; 2570 mat_ilen = c->ilen + i; 2571 for (k=kstart; k<kend; k++) { 2572 if ((tcol=smap[a->j[k]])) { 2573 *mat_j++ = tcol - 1; 2574 *mat_a++ = a->a[k]; 2575 (*mat_ilen)++; 2576 2577 } 2578 } 2579 } 2580 /* Free work space */ 2581 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2582 ierr = PetscFree(smap);CHKERRQ(ierr); 2583 ierr = PetscFree(lens);CHKERRQ(ierr); 2584 } 2585 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2586 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2587 2588 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2589 *B = C; 2590 PetscFunctionReturn(0); 2591 } 2592 2593 #undef __FUNCT__ 2594 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2595 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2596 { 2597 PetscErrorCode ierr; 2598 Mat B; 2599 2600 PetscFunctionBegin; 2601 if (scall == MAT_INITIAL_MATRIX) { 2602 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2603 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2604 ierr = MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs);CHKERRQ(ierr); 2605 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2606 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2607 *subMat = B; 2608 } else { 2609 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2610 } 2611 PetscFunctionReturn(0); 2612 } 2613 2614 #undef __FUNCT__ 2615 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2616 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2617 { 2618 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2619 PetscErrorCode ierr; 2620 Mat outA; 2621 PetscBool row_identity,col_identity; 2622 2623 PetscFunctionBegin; 2624 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2625 2626 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2627 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2628 2629 outA = inA; 2630 outA->factortype = MAT_FACTOR_LU; 2631 2632 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2633 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2634 2635 a->row = row; 2636 2637 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2638 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2639 2640 a->col = col; 2641 2642 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2643 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2644 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2645 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2646 2647 if (!a->solve_work) { /* this matrix may have been factored before */ 2648 ierr = PetscMalloc1((inA->rmap->n+1),&a->solve_work);CHKERRQ(ierr); 2649 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2650 } 2651 2652 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2653 if (row_identity && col_identity) { 2654 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2655 } else { 2656 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2657 } 2658 PetscFunctionReturn(0); 2659 } 2660 2661 #undef __FUNCT__ 2662 #define __FUNCT__ "MatScale_SeqAIJ" 2663 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2664 { 2665 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2666 PetscScalar oalpha = alpha; 2667 PetscErrorCode ierr; 2668 PetscBLASInt one = 1,bnz; 2669 2670 PetscFunctionBegin; 2671 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2672 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2673 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2674 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2675 PetscFunctionReturn(0); 2676 } 2677 2678 #undef __FUNCT__ 2679 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2680 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2681 { 2682 PetscErrorCode ierr; 2683 PetscInt i; 2684 2685 PetscFunctionBegin; 2686 if (scall == MAT_INITIAL_MATRIX) { 2687 ierr = PetscMalloc1((n+1),B);CHKERRQ(ierr); 2688 } 2689 2690 for (i=0; i<n; i++) { 2691 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2692 } 2693 PetscFunctionReturn(0); 2694 } 2695 2696 #undef __FUNCT__ 2697 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2698 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2699 { 2700 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2701 PetscErrorCode ierr; 2702 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2703 const PetscInt *idx; 2704 PetscInt start,end,*ai,*aj; 2705 PetscBT table; 2706 2707 PetscFunctionBegin; 2708 m = A->rmap->n; 2709 ai = a->i; 2710 aj = a->j; 2711 2712 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2713 2714 ierr = PetscMalloc1((m+1),&nidx);CHKERRQ(ierr); 2715 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2716 2717 for (i=0; i<is_max; i++) { 2718 /* Initialize the two local arrays */ 2719 isz = 0; 2720 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2721 2722 /* Extract the indices, assume there can be duplicate entries */ 2723 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2724 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2725 2726 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2727 for (j=0; j<n; ++j) { 2728 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2729 } 2730 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2731 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2732 2733 k = 0; 2734 for (j=0; j<ov; j++) { /* for each overlap */ 2735 n = isz; 2736 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2737 row = nidx[k]; 2738 start = ai[row]; 2739 end = ai[row+1]; 2740 for (l = start; l<end; l++) { 2741 val = aj[l]; 2742 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2743 } 2744 } 2745 } 2746 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2747 } 2748 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2749 ierr = PetscFree(nidx);CHKERRQ(ierr); 2750 PetscFunctionReturn(0); 2751 } 2752 2753 /* -------------------------------------------------------------- */ 2754 #undef __FUNCT__ 2755 #define __FUNCT__ "MatPermute_SeqAIJ" 2756 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2757 { 2758 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2759 PetscErrorCode ierr; 2760 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2761 const PetscInt *row,*col; 2762 PetscInt *cnew,j,*lens; 2763 IS icolp,irowp; 2764 PetscInt *cwork = NULL; 2765 PetscScalar *vwork = NULL; 2766 2767 PetscFunctionBegin; 2768 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2769 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2770 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2771 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2772 2773 /* determine lengths of permuted rows */ 2774 ierr = PetscMalloc1((m+1),&lens);CHKERRQ(ierr); 2775 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2776 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2777 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2778 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2779 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2780 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2781 ierr = PetscFree(lens);CHKERRQ(ierr); 2782 2783 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2784 for (i=0; i<m; i++) { 2785 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2786 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2787 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2788 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2789 } 2790 ierr = PetscFree(cnew);CHKERRQ(ierr); 2791 2792 (*B)->assembled = PETSC_FALSE; 2793 2794 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2795 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2796 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2797 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2798 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2799 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2800 PetscFunctionReturn(0); 2801 } 2802 2803 #undef __FUNCT__ 2804 #define __FUNCT__ "MatCopy_SeqAIJ" 2805 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2806 { 2807 PetscErrorCode ierr; 2808 2809 PetscFunctionBegin; 2810 /* If the two matrices have the same copy implementation, use fast copy. */ 2811 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2812 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2813 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2814 2815 if (a->i[A->rmap->n] != b->i[B->rmap->n]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2816 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2817 } else { 2818 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2819 } 2820 PetscFunctionReturn(0); 2821 } 2822 2823 #undef __FUNCT__ 2824 #define __FUNCT__ "MatSetUp_SeqAIJ" 2825 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2826 { 2827 PetscErrorCode ierr; 2828 2829 PetscFunctionBegin; 2830 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2831 PetscFunctionReturn(0); 2832 } 2833 2834 #undef __FUNCT__ 2835 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2836 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2837 { 2838 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2839 2840 PetscFunctionBegin; 2841 *array = a->a; 2842 PetscFunctionReturn(0); 2843 } 2844 2845 #undef __FUNCT__ 2846 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2847 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2848 { 2849 PetscFunctionBegin; 2850 PetscFunctionReturn(0); 2851 } 2852 2853 /* 2854 Computes the number of nonzeros per row needed for preallocation when X and Y 2855 have different nonzero structure. 2856 */ 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2859 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2860 { 2861 PetscInt i,m=Y->rmap->N; 2862 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2863 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2864 const PetscInt *xi = x->i,*yi = y->i; 2865 2866 PetscFunctionBegin; 2867 /* Set the number of nonzeros in the new matrix */ 2868 for (i=0; i<m; i++) { 2869 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2870 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2871 nnz[i] = 0; 2872 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2873 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2874 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2875 nnz[i]++; 2876 } 2877 for (; k<nzy; k++) nnz[i]++; 2878 } 2879 PetscFunctionReturn(0); 2880 } 2881 2882 #undef __FUNCT__ 2883 #define __FUNCT__ "MatAXPY_SeqAIJ" 2884 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2885 { 2886 PetscErrorCode ierr; 2887 PetscInt i; 2888 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2889 PetscBLASInt one=1,bnz; 2890 2891 PetscFunctionBegin; 2892 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2893 if (str == SAME_NONZERO_PATTERN) { 2894 PetscScalar alpha = a; 2895 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2896 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2897 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2898 if (y->xtoy && y->XtoY != X) { 2899 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2900 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2901 } 2902 if (!y->xtoy) { /* get xtoy */ 2903 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,NULL, y->i,y->j,NULL, &y->xtoy);CHKERRQ(ierr); 2904 y->XtoY = X; 2905 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2906 } 2907 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2908 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(double)(PetscReal)(x->nz)/(y->nz+1));CHKERRQ(ierr); 2909 } else { 2910 Mat B; 2911 PetscInt *nnz; 2912 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2913 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2914 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2915 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2916 ierr = MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);CHKERRQ(ierr); 2917 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2918 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2919 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2920 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2921 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2922 ierr = PetscFree(nnz);CHKERRQ(ierr); 2923 } 2924 PetscFunctionReturn(0); 2925 } 2926 2927 #undef __FUNCT__ 2928 #define __FUNCT__ "MatConjugate_SeqAIJ" 2929 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2930 { 2931 #if defined(PETSC_USE_COMPLEX) 2932 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2933 PetscInt i,nz; 2934 PetscScalar *a; 2935 2936 PetscFunctionBegin; 2937 nz = aij->nz; 2938 a = aij->a; 2939 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2940 #else 2941 PetscFunctionBegin; 2942 #endif 2943 PetscFunctionReturn(0); 2944 } 2945 2946 #undef __FUNCT__ 2947 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2948 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2949 { 2950 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2951 PetscErrorCode ierr; 2952 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2953 PetscReal atmp; 2954 PetscScalar *x; 2955 MatScalar *aa; 2956 2957 PetscFunctionBegin; 2958 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2959 aa = a->a; 2960 ai = a->i; 2961 aj = a->j; 2962 2963 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2964 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2965 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2966 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2967 for (i=0; i<m; i++) { 2968 ncols = ai[1] - ai[0]; ai++; 2969 x[i] = 0.0; 2970 for (j=0; j<ncols; j++) { 2971 atmp = PetscAbsScalar(*aa); 2972 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2973 aa++; aj++; 2974 } 2975 } 2976 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2977 PetscFunctionReturn(0); 2978 } 2979 2980 #undef __FUNCT__ 2981 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2982 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2983 { 2984 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2985 PetscErrorCode ierr; 2986 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2987 PetscScalar *x; 2988 MatScalar *aa; 2989 2990 PetscFunctionBegin; 2991 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2992 aa = a->a; 2993 ai = a->i; 2994 aj = a->j; 2995 2996 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2997 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2998 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2999 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3000 for (i=0; i<m; i++) { 3001 ncols = ai[1] - ai[0]; ai++; 3002 if (ncols == A->cmap->n) { /* row is dense */ 3003 x[i] = *aa; if (idx) idx[i] = 0; 3004 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3005 x[i] = 0.0; 3006 if (idx) { 3007 idx[i] = 0; /* in case ncols is zero */ 3008 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 3009 if (aj[j] > j) { 3010 idx[i] = j; 3011 break; 3012 } 3013 } 3014 } 3015 } 3016 for (j=0; j<ncols; j++) { 3017 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3018 aa++; aj++; 3019 } 3020 } 3021 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3022 PetscFunctionReturn(0); 3023 } 3024 3025 #undef __FUNCT__ 3026 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 3027 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3028 { 3029 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3030 PetscErrorCode ierr; 3031 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3032 PetscReal atmp; 3033 PetscScalar *x; 3034 MatScalar *aa; 3035 3036 PetscFunctionBegin; 3037 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3038 aa = a->a; 3039 ai = a->i; 3040 aj = a->j; 3041 3042 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3043 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3044 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3045 if (n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %d vs. %d rows", A->rmap->n, n); 3046 for (i=0; i<m; i++) { 3047 ncols = ai[1] - ai[0]; ai++; 3048 if (ncols) { 3049 /* Get first nonzero */ 3050 for (j = 0; j < ncols; j++) { 3051 atmp = PetscAbsScalar(aa[j]); 3052 if (atmp > 1.0e-12) { 3053 x[i] = atmp; 3054 if (idx) idx[i] = aj[j]; 3055 break; 3056 } 3057 } 3058 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3059 } else { 3060 x[i] = 0.0; if (idx) idx[i] = 0; 3061 } 3062 for (j = 0; j < ncols; j++) { 3063 atmp = PetscAbsScalar(*aa); 3064 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3065 aa++; aj++; 3066 } 3067 } 3068 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3069 PetscFunctionReturn(0); 3070 } 3071 3072 #undef __FUNCT__ 3073 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 3074 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3075 { 3076 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3077 PetscErrorCode ierr; 3078 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3079 PetscScalar *x; 3080 MatScalar *aa; 3081 3082 PetscFunctionBegin; 3083 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3084 aa = a->a; 3085 ai = a->i; 3086 aj = a->j; 3087 3088 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3089 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3090 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3091 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3092 for (i=0; i<m; i++) { 3093 ncols = ai[1] - ai[0]; ai++; 3094 if (ncols == A->cmap->n) { /* row is dense */ 3095 x[i] = *aa; if (idx) idx[i] = 0; 3096 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3097 x[i] = 0.0; 3098 if (idx) { /* find first implicit 0.0 in the row */ 3099 idx[i] = 0; /* in case ncols is zero */ 3100 for (j=0; j<ncols; j++) { 3101 if (aj[j] > j) { 3102 idx[i] = j; 3103 break; 3104 } 3105 } 3106 } 3107 } 3108 for (j=0; j<ncols; j++) { 3109 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3110 aa++; aj++; 3111 } 3112 } 3113 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3114 PetscFunctionReturn(0); 3115 } 3116 3117 #include <petscblaslapack.h> 3118 #include <petsc-private/kernels/blockinvert.h> 3119 3120 #undef __FUNCT__ 3121 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 3122 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3123 { 3124 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3125 PetscErrorCode ierr; 3126 PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3127 MatScalar *diag,work[25],*v_work; 3128 PetscReal shift = 0.0; 3129 3130 PetscFunctionBegin; 3131 if (a->ibdiagvalid) { 3132 if (values) *values = a->ibdiag; 3133 PetscFunctionReturn(0); 3134 } 3135 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3136 if (!a->ibdiag) { 3137 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3138 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3139 } 3140 diag = a->ibdiag; 3141 if (values) *values = a->ibdiag; 3142 /* factor and invert each block */ 3143 switch (bs) { 3144 case 1: 3145 for (i=0; i<mbs; i++) { 3146 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3147 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3148 } 3149 break; 3150 case 2: 3151 for (i=0; i<mbs; i++) { 3152 ij[0] = 2*i; ij[1] = 2*i + 1; 3153 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3154 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 3155 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3156 diag += 4; 3157 } 3158 break; 3159 case 3: 3160 for (i=0; i<mbs; i++) { 3161 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3162 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3163 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 3164 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3165 diag += 9; 3166 } 3167 break; 3168 case 4: 3169 for (i=0; i<mbs; i++) { 3170 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3171 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3172 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 3173 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3174 diag += 16; 3175 } 3176 break; 3177 case 5: 3178 for (i=0; i<mbs; i++) { 3179 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3180 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3181 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 3182 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3183 diag += 25; 3184 } 3185 break; 3186 case 6: 3187 for (i=0; i<mbs; i++) { 3188 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3189 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3190 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 3191 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3192 diag += 36; 3193 } 3194 break; 3195 case 7: 3196 for (i=0; i<mbs; i++) { 3197 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3198 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3199 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3200 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3201 diag += 49; 3202 } 3203 break; 3204 default: 3205 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3206 for (i=0; i<mbs; i++) { 3207 for (j=0; j<bs; j++) { 3208 IJ[j] = bs*i + j; 3209 } 3210 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3211 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3212 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3213 diag += bs2; 3214 } 3215 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3216 } 3217 a->ibdiagvalid = PETSC_TRUE; 3218 PetscFunctionReturn(0); 3219 } 3220 3221 #undef __FUNCT__ 3222 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3223 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3224 { 3225 PetscErrorCode ierr; 3226 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3227 PetscScalar a; 3228 PetscInt m,n,i,j,col; 3229 3230 PetscFunctionBegin; 3231 if (!x->assembled) { 3232 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3233 for (i=0; i<m; i++) { 3234 for (j=0; j<aij->imax[i]; j++) { 3235 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3236 col = (PetscInt)(n*PetscRealPart(a)); 3237 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3238 } 3239 } 3240 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3241 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3242 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3243 PetscFunctionReturn(0); 3244 } 3245 3246 /* -------------------------------------------------------------------*/ 3247 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3248 MatGetRow_SeqAIJ, 3249 MatRestoreRow_SeqAIJ, 3250 MatMult_SeqAIJ, 3251 /* 4*/ MatMultAdd_SeqAIJ, 3252 MatMultTranspose_SeqAIJ, 3253 MatMultTransposeAdd_SeqAIJ, 3254 0, 3255 0, 3256 0, 3257 /* 10*/ 0, 3258 MatLUFactor_SeqAIJ, 3259 0, 3260 MatSOR_SeqAIJ, 3261 MatTranspose_SeqAIJ, 3262 /*1 5*/ MatGetInfo_SeqAIJ, 3263 MatEqual_SeqAIJ, 3264 MatGetDiagonal_SeqAIJ, 3265 MatDiagonalScale_SeqAIJ, 3266 MatNorm_SeqAIJ, 3267 /* 20*/ 0, 3268 MatAssemblyEnd_SeqAIJ, 3269 MatSetOption_SeqAIJ, 3270 MatZeroEntries_SeqAIJ, 3271 /* 24*/ MatZeroRows_SeqAIJ, 3272 0, 3273 0, 3274 0, 3275 0, 3276 /* 29*/ MatSetUp_SeqAIJ, 3277 0, 3278 0, 3279 0, 3280 0, 3281 /* 34*/ MatDuplicate_SeqAIJ, 3282 0, 3283 0, 3284 MatILUFactor_SeqAIJ, 3285 0, 3286 /* 39*/ MatAXPY_SeqAIJ, 3287 MatGetSubMatrices_SeqAIJ, 3288 MatIncreaseOverlap_SeqAIJ, 3289 MatGetValues_SeqAIJ, 3290 MatCopy_SeqAIJ, 3291 /* 44*/ MatGetRowMax_SeqAIJ, 3292 MatScale_SeqAIJ, 3293 0, 3294 MatDiagonalSet_SeqAIJ, 3295 MatZeroRowsColumns_SeqAIJ, 3296 /* 49*/ MatSetRandom_SeqAIJ, 3297 MatGetRowIJ_SeqAIJ, 3298 MatRestoreRowIJ_SeqAIJ, 3299 MatGetColumnIJ_SeqAIJ, 3300 MatRestoreColumnIJ_SeqAIJ, 3301 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3302 0, 3303 0, 3304 MatPermute_SeqAIJ, 3305 0, 3306 /* 59*/ 0, 3307 MatDestroy_SeqAIJ, 3308 MatView_SeqAIJ, 3309 0, 3310 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3311 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3312 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3313 0, 3314 0, 3315 0, 3316 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3317 MatGetRowMinAbs_SeqAIJ, 3318 0, 3319 MatSetColoring_SeqAIJ, 3320 0, 3321 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3322 MatFDColoringApply_AIJ, 3323 0, 3324 0, 3325 0, 3326 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3327 0, 3328 0, 3329 0, 3330 MatLoad_SeqAIJ, 3331 /* 84*/ MatIsSymmetric_SeqAIJ, 3332 MatIsHermitian_SeqAIJ, 3333 0, 3334 0, 3335 0, 3336 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3337 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3338 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3339 MatPtAP_SeqAIJ_SeqAIJ, 3340 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3341 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3342 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3343 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3344 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3345 0, 3346 /* 99*/ 0, 3347 0, 3348 0, 3349 MatConjugate_SeqAIJ, 3350 0, 3351 /*104*/ MatSetValuesRow_SeqAIJ, 3352 MatRealPart_SeqAIJ, 3353 MatImaginaryPart_SeqAIJ, 3354 0, 3355 0, 3356 /*109*/ MatMatSolve_SeqAIJ, 3357 0, 3358 MatGetRowMin_SeqAIJ, 3359 0, 3360 MatMissingDiagonal_SeqAIJ, 3361 /*114*/ 0, 3362 0, 3363 0, 3364 0, 3365 0, 3366 /*119*/ 0, 3367 0, 3368 0, 3369 0, 3370 MatGetMultiProcBlock_SeqAIJ, 3371 /*124*/ MatFindNonzeroRows_SeqAIJ, 3372 MatGetColumnNorms_SeqAIJ, 3373 MatInvertBlockDiagonal_SeqAIJ, 3374 0, 3375 0, 3376 /*129*/ 0, 3377 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3378 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3379 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3380 MatTransposeColoringCreate_SeqAIJ, 3381 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3382 MatTransColoringApplyDenToSp_SeqAIJ, 3383 MatRARt_SeqAIJ_SeqAIJ, 3384 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3385 MatRARtNumeric_SeqAIJ_SeqAIJ, 3386 /*139*/0, 3387 0, 3388 0, 3389 MatFDColoringSetUp_SeqXAIJ, 3390 MatFindOffBlockDiagonalEntries_SeqAIJ 3391 }; 3392 3393 #undef __FUNCT__ 3394 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3395 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3396 { 3397 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3398 PetscInt i,nz,n; 3399 3400 PetscFunctionBegin; 3401 nz = aij->maxnz; 3402 n = mat->rmap->n; 3403 for (i=0; i<nz; i++) { 3404 aij->j[i] = indices[i]; 3405 } 3406 aij->nz = nz; 3407 for (i=0; i<n; i++) { 3408 aij->ilen[i] = aij->imax[i]; 3409 } 3410 PetscFunctionReturn(0); 3411 } 3412 3413 #undef __FUNCT__ 3414 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3415 /*@ 3416 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3417 in the matrix. 3418 3419 Input Parameters: 3420 + mat - the SeqAIJ matrix 3421 - indices - the column indices 3422 3423 Level: advanced 3424 3425 Notes: 3426 This can be called if you have precomputed the nonzero structure of the 3427 matrix and want to provide it to the matrix object to improve the performance 3428 of the MatSetValues() operation. 3429 3430 You MUST have set the correct numbers of nonzeros per row in the call to 3431 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3432 3433 MUST be called before any calls to MatSetValues(); 3434 3435 The indices should start with zero, not one. 3436 3437 @*/ 3438 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3439 { 3440 PetscErrorCode ierr; 3441 3442 PetscFunctionBegin; 3443 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3444 PetscValidPointer(indices,2); 3445 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3446 PetscFunctionReturn(0); 3447 } 3448 3449 /* ----------------------------------------------------------------------------------------*/ 3450 3451 #undef __FUNCT__ 3452 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3453 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3454 { 3455 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3456 PetscErrorCode ierr; 3457 size_t nz = aij->i[mat->rmap->n]; 3458 3459 PetscFunctionBegin; 3460 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3461 3462 /* allocate space for values if not already there */ 3463 if (!aij->saved_values) { 3464 ierr = PetscMalloc1((nz+1),&aij->saved_values);CHKERRQ(ierr); 3465 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3466 } 3467 3468 /* copy values over */ 3469 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3470 PetscFunctionReturn(0); 3471 } 3472 3473 #undef __FUNCT__ 3474 #define __FUNCT__ "MatStoreValues" 3475 /*@ 3476 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3477 example, reuse of the linear part of a Jacobian, while recomputing the 3478 nonlinear portion. 3479 3480 Collect on Mat 3481 3482 Input Parameters: 3483 . mat - the matrix (currently only AIJ matrices support this option) 3484 3485 Level: advanced 3486 3487 Common Usage, with SNESSolve(): 3488 $ Create Jacobian matrix 3489 $ Set linear terms into matrix 3490 $ Apply boundary conditions to matrix, at this time matrix must have 3491 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3492 $ boundary conditions again will not change the nonzero structure 3493 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3494 $ ierr = MatStoreValues(mat); 3495 $ Call SNESSetJacobian() with matrix 3496 $ In your Jacobian routine 3497 $ ierr = MatRetrieveValues(mat); 3498 $ Set nonlinear terms in matrix 3499 3500 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3501 $ // build linear portion of Jacobian 3502 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3503 $ ierr = MatStoreValues(mat); 3504 $ loop over nonlinear iterations 3505 $ ierr = MatRetrieveValues(mat); 3506 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3507 $ // call MatAssemblyBegin/End() on matrix 3508 $ Solve linear system with Jacobian 3509 $ endloop 3510 3511 Notes: 3512 Matrix must already be assemblied before calling this routine 3513 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3514 calling this routine. 3515 3516 When this is called multiple times it overwrites the previous set of stored values 3517 and does not allocated additional space. 3518 3519 .seealso: MatRetrieveValues() 3520 3521 @*/ 3522 PetscErrorCode MatStoreValues(Mat mat) 3523 { 3524 PetscErrorCode ierr; 3525 3526 PetscFunctionBegin; 3527 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3528 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3529 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3530 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3531 PetscFunctionReturn(0); 3532 } 3533 3534 #undef __FUNCT__ 3535 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3536 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3537 { 3538 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3539 PetscErrorCode ierr; 3540 PetscInt nz = aij->i[mat->rmap->n]; 3541 3542 PetscFunctionBegin; 3543 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3544 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3545 /* copy values over */ 3546 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3547 PetscFunctionReturn(0); 3548 } 3549 3550 #undef __FUNCT__ 3551 #define __FUNCT__ "MatRetrieveValues" 3552 /*@ 3553 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3554 example, reuse of the linear part of a Jacobian, while recomputing the 3555 nonlinear portion. 3556 3557 Collect on Mat 3558 3559 Input Parameters: 3560 . mat - the matrix (currently on AIJ matrices support this option) 3561 3562 Level: advanced 3563 3564 .seealso: MatStoreValues() 3565 3566 @*/ 3567 PetscErrorCode MatRetrieveValues(Mat mat) 3568 { 3569 PetscErrorCode ierr; 3570 3571 PetscFunctionBegin; 3572 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3573 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3574 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3575 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3576 PetscFunctionReturn(0); 3577 } 3578 3579 3580 /* --------------------------------------------------------------------------------*/ 3581 #undef __FUNCT__ 3582 #define __FUNCT__ "MatCreateSeqAIJ" 3583 /*@C 3584 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3585 (the default parallel PETSc format). For good matrix assembly performance 3586 the user should preallocate the matrix storage by setting the parameter nz 3587 (or the array nnz). By setting these parameters accurately, performance 3588 during matrix assembly can be increased by more than a factor of 50. 3589 3590 Collective on MPI_Comm 3591 3592 Input Parameters: 3593 + comm - MPI communicator, set to PETSC_COMM_SELF 3594 . m - number of rows 3595 . n - number of columns 3596 . nz - number of nonzeros per row (same for all rows) 3597 - nnz - array containing the number of nonzeros in the various rows 3598 (possibly different for each row) or NULL 3599 3600 Output Parameter: 3601 . A - the matrix 3602 3603 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3604 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3605 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3606 3607 Notes: 3608 If nnz is given then nz is ignored 3609 3610 The AIJ format (also called the Yale sparse matrix format or 3611 compressed row storage), is fully compatible with standard Fortran 77 3612 storage. That is, the stored row and column indices can begin at 3613 either one (as in Fortran) or zero. See the users' manual for details. 3614 3615 Specify the preallocated storage with either nz or nnz (not both). 3616 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3617 allocation. For large problems you MUST preallocate memory or you 3618 will get TERRIBLE performance, see the users' manual chapter on matrices. 3619 3620 By default, this format uses inodes (identical nodes) when possible, to 3621 improve numerical efficiency of matrix-vector products and solves. We 3622 search for consecutive rows with the same nonzero structure, thereby 3623 reusing matrix information to achieve increased efficiency. 3624 3625 Options Database Keys: 3626 + -mat_no_inode - Do not use inodes 3627 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3628 3629 Level: intermediate 3630 3631 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3632 3633 @*/ 3634 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3635 { 3636 PetscErrorCode ierr; 3637 3638 PetscFunctionBegin; 3639 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3640 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3641 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3642 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3643 PetscFunctionReturn(0); 3644 } 3645 3646 #undef __FUNCT__ 3647 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3648 /*@C 3649 MatSeqAIJSetPreallocation - For good matrix assembly performance 3650 the user should preallocate the matrix storage by setting the parameter nz 3651 (or the array nnz). By setting these parameters accurately, performance 3652 during matrix assembly can be increased by more than a factor of 50. 3653 3654 Collective on MPI_Comm 3655 3656 Input Parameters: 3657 + B - The matrix-free 3658 . nz - number of nonzeros per row (same for all rows) 3659 - nnz - array containing the number of nonzeros in the various rows 3660 (possibly different for each row) or NULL 3661 3662 Notes: 3663 If nnz is given then nz is ignored 3664 3665 The AIJ format (also called the Yale sparse matrix format or 3666 compressed row storage), is fully compatible with standard Fortran 77 3667 storage. That is, the stored row and column indices can begin at 3668 either one (as in Fortran) or zero. See the users' manual for details. 3669 3670 Specify the preallocated storage with either nz or nnz (not both). 3671 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3672 allocation. For large problems you MUST preallocate memory or you 3673 will get TERRIBLE performance, see the users' manual chapter on matrices. 3674 3675 You can call MatGetInfo() to get information on how effective the preallocation was; 3676 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3677 You can also run with the option -info and look for messages with the string 3678 malloc in them to see if additional memory allocation was needed. 3679 3680 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3681 entries or columns indices 3682 3683 By default, this format uses inodes (identical nodes) when possible, to 3684 improve numerical efficiency of matrix-vector products and solves. We 3685 search for consecutive rows with the same nonzero structure, thereby 3686 reusing matrix information to achieve increased efficiency. 3687 3688 Options Database Keys: 3689 + -mat_no_inode - Do not use inodes 3690 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3691 - -mat_aij_oneindex - Internally use indexing starting at 1 3692 rather than 0. Note that when calling MatSetValues(), 3693 the user still MUST index entries starting at 0! 3694 3695 Level: intermediate 3696 3697 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3698 3699 @*/ 3700 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3701 { 3702 PetscErrorCode ierr; 3703 3704 PetscFunctionBegin; 3705 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3706 PetscValidType(B,1); 3707 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3708 PetscFunctionReturn(0); 3709 } 3710 3711 #undef __FUNCT__ 3712 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3713 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3714 { 3715 Mat_SeqAIJ *b; 3716 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3717 PetscErrorCode ierr; 3718 PetscInt i; 3719 3720 PetscFunctionBegin; 3721 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3722 if (nz == MAT_SKIP_ALLOCATION) { 3723 skipallocation = PETSC_TRUE; 3724 nz = 0; 3725 } 3726 3727 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3728 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3729 3730 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3731 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3732 if (nnz) { 3733 for (i=0; i<B->rmap->n; i++) { 3734 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 3735 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3736 } 3737 } 3738 3739 B->preallocated = PETSC_TRUE; 3740 3741 b = (Mat_SeqAIJ*)B->data; 3742 3743 if (!skipallocation) { 3744 if (!b->imax) { 3745 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3746 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3747 } 3748 if (!nnz) { 3749 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3750 else if (nz < 0) nz = 1; 3751 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3752 nz = nz*B->rmap->n; 3753 } else { 3754 nz = 0; 3755 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3756 } 3757 /* b->ilen will count nonzeros in each row so far. */ 3758 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3759 3760 /* allocate the matrix space */ 3761 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3762 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3763 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3764 b->i[0] = 0; 3765 for (i=1; i<B->rmap->n+1; i++) { 3766 b->i[i] = b->i[i-1] + b->imax[i-1]; 3767 } 3768 b->singlemalloc = PETSC_TRUE; 3769 b->free_a = PETSC_TRUE; 3770 b->free_ij = PETSC_TRUE; 3771 #if defined(PETSC_THREADCOMM_ACTIVE) 3772 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3773 #endif 3774 } else { 3775 b->free_a = PETSC_FALSE; 3776 b->free_ij = PETSC_FALSE; 3777 } 3778 3779 b->nz = 0; 3780 b->maxnz = nz; 3781 B->info.nz_unneeded = (double)b->maxnz; 3782 if (realalloc) { 3783 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3784 } 3785 PetscFunctionReturn(0); 3786 } 3787 3788 #undef __FUNCT__ 3789 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3790 /*@ 3791 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3792 3793 Input Parameters: 3794 + B - the matrix 3795 . i - the indices into j for the start of each row (starts with zero) 3796 . j - the column indices for each row (starts with zero) these must be sorted for each row 3797 - v - optional values in the matrix 3798 3799 Level: developer 3800 3801 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3802 3803 .keywords: matrix, aij, compressed row, sparse, sequential 3804 3805 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3806 @*/ 3807 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3808 { 3809 PetscErrorCode ierr; 3810 3811 PetscFunctionBegin; 3812 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3813 PetscValidType(B,1); 3814 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3815 PetscFunctionReturn(0); 3816 } 3817 3818 #undef __FUNCT__ 3819 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3820 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3821 { 3822 PetscInt i; 3823 PetscInt m,n; 3824 PetscInt nz; 3825 PetscInt *nnz, nz_max = 0; 3826 PetscScalar *values; 3827 PetscErrorCode ierr; 3828 3829 PetscFunctionBegin; 3830 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3831 3832 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3833 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3834 3835 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3836 ierr = PetscMalloc1((m+1), &nnz);CHKERRQ(ierr); 3837 for (i = 0; i < m; i++) { 3838 nz = Ii[i+1]- Ii[i]; 3839 nz_max = PetscMax(nz_max, nz); 3840 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3841 nnz[i] = nz; 3842 } 3843 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3844 ierr = PetscFree(nnz);CHKERRQ(ierr); 3845 3846 if (v) { 3847 values = (PetscScalar*) v; 3848 } else { 3849 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3850 } 3851 3852 for (i = 0; i < m; i++) { 3853 nz = Ii[i+1] - Ii[i]; 3854 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3855 } 3856 3857 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3858 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3859 3860 if (!v) { 3861 ierr = PetscFree(values);CHKERRQ(ierr); 3862 } 3863 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3864 PetscFunctionReturn(0); 3865 } 3866 3867 #include <../src/mat/impls/dense/seq/dense.h> 3868 #include <petsc-private/kernels/petscaxpy.h> 3869 3870 #undef __FUNCT__ 3871 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3872 /* 3873 Computes (B'*A')' since computing B*A directly is untenable 3874 3875 n p p 3876 ( ) ( ) ( ) 3877 m ( A ) * n ( B ) = m ( C ) 3878 ( ) ( ) ( ) 3879 3880 */ 3881 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3882 { 3883 PetscErrorCode ierr; 3884 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3885 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3886 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3887 PetscInt i,n,m,q,p; 3888 const PetscInt *ii,*idx; 3889 const PetscScalar *b,*a,*a_q; 3890 PetscScalar *c,*c_q; 3891 3892 PetscFunctionBegin; 3893 m = A->rmap->n; 3894 n = A->cmap->n; 3895 p = B->cmap->n; 3896 a = sub_a->v; 3897 b = sub_b->a; 3898 c = sub_c->v; 3899 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3900 3901 ii = sub_b->i; 3902 idx = sub_b->j; 3903 for (i=0; i<n; i++) { 3904 q = ii[i+1] - ii[i]; 3905 while (q-->0) { 3906 c_q = c + m*(*idx); 3907 a_q = a + m*i; 3908 PetscKernelAXPY(c_q,*b,a_q,m); 3909 idx++; 3910 b++; 3911 } 3912 } 3913 PetscFunctionReturn(0); 3914 } 3915 3916 #undef __FUNCT__ 3917 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3918 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3919 { 3920 PetscErrorCode ierr; 3921 PetscInt m=A->rmap->n,n=B->cmap->n; 3922 Mat Cmat; 3923 3924 PetscFunctionBegin; 3925 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3926 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3927 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3928 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3929 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3930 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3931 3932 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3933 3934 *C = Cmat; 3935 PetscFunctionReturn(0); 3936 } 3937 3938 /* ----------------------------------------------------------------*/ 3939 #undef __FUNCT__ 3940 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3941 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3942 { 3943 PetscErrorCode ierr; 3944 3945 PetscFunctionBegin; 3946 if (scall == MAT_INITIAL_MATRIX) { 3947 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3948 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3949 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3950 } 3951 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3952 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3953 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3954 PetscFunctionReturn(0); 3955 } 3956 3957 3958 /*MC 3959 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3960 based on compressed sparse row format. 3961 3962 Options Database Keys: 3963 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3964 3965 Level: beginner 3966 3967 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3968 M*/ 3969 3970 /*MC 3971 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3972 3973 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3974 and MATMPIAIJ otherwise. As a result, for single process communicators, 3975 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3976 for communicators controlling multiple processes. It is recommended that you call both of 3977 the above preallocation routines for simplicity. 3978 3979 Options Database Keys: 3980 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3981 3982 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3983 enough exist. 3984 3985 Level: beginner 3986 3987 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3988 M*/ 3989 3990 /*MC 3991 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3992 3993 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3994 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3995 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3996 for communicators controlling multiple processes. It is recommended that you call both of 3997 the above preallocation routines for simplicity. 3998 3999 Options Database Keys: 4000 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4001 4002 Level: beginner 4003 4004 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4005 M*/ 4006 4007 #if defined(PETSC_HAVE_PASTIX) 4008 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 4009 #endif 4010 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 4011 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*); 4012 #endif 4013 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4014 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 4015 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 4016 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*); 4017 #if defined(PETSC_HAVE_MUMPS) 4018 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 4019 #endif 4020 #if defined(PETSC_HAVE_SUPERLU) 4021 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 4022 #endif 4023 #if defined(PETSC_HAVE_SUPERLU_DIST) 4024 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 4025 #endif 4026 #if defined(PETSC_HAVE_UMFPACK) 4027 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 4028 #endif 4029 #if defined(PETSC_HAVE_CHOLMOD) 4030 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 4031 #endif 4032 #if defined(PETSC_HAVE_LUSOL) 4033 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 4034 #endif 4035 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4036 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 4037 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 4038 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 4039 #endif 4040 #if defined(PETSC_HAVE_CLIQUE) 4041 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 4042 #endif 4043 4044 4045 #undef __FUNCT__ 4046 #define __FUNCT__ "MatSeqAIJGetArray" 4047 /*@C 4048 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 4049 4050 Not Collective 4051 4052 Input Parameter: 4053 . mat - a MATSEQDENSE matrix 4054 4055 Output Parameter: 4056 . array - pointer to the data 4057 4058 Level: intermediate 4059 4060 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4061 @*/ 4062 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4063 { 4064 PetscErrorCode ierr; 4065 4066 PetscFunctionBegin; 4067 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4068 PetscFunctionReturn(0); 4069 } 4070 4071 #undef __FUNCT__ 4072 #define __FUNCT__ "MatSeqAIJRestoreArray" 4073 /*@C 4074 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 4075 4076 Not Collective 4077 4078 Input Parameters: 4079 . mat - a MATSEQDENSE matrix 4080 . array - pointer to the data 4081 4082 Level: intermediate 4083 4084 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4085 @*/ 4086 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4087 { 4088 PetscErrorCode ierr; 4089 4090 PetscFunctionBegin; 4091 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4092 PetscFunctionReturn(0); 4093 } 4094 4095 #undef __FUNCT__ 4096 #define __FUNCT__ "MatCreate_SeqAIJ" 4097 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4098 { 4099 Mat_SeqAIJ *b; 4100 PetscErrorCode ierr; 4101 PetscMPIInt size; 4102 4103 PetscFunctionBegin; 4104 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4105 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4106 4107 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4108 4109 B->data = (void*)b; 4110 4111 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4112 4113 b->row = 0; 4114 b->col = 0; 4115 b->icol = 0; 4116 b->reallocs = 0; 4117 b->ignorezeroentries = PETSC_FALSE; 4118 b->roworiented = PETSC_TRUE; 4119 b->nonew = 0; 4120 b->diag = 0; 4121 b->solve_work = 0; 4122 B->spptr = 0; 4123 b->saved_values = 0; 4124 b->idiag = 0; 4125 b->mdiag = 0; 4126 b->ssor_work = 0; 4127 b->omega = 1.0; 4128 b->fshift = 0.0; 4129 b->idiagvalid = PETSC_FALSE; 4130 b->ibdiagvalid = PETSC_FALSE; 4131 b->keepnonzeropattern = PETSC_FALSE; 4132 b->xtoy = 0; 4133 b->XtoY = 0; 4134 B->same_nonzero = PETSC_FALSE; 4135 4136 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4137 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4138 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4139 4140 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4141 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 4142 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4143 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4144 #endif 4145 #if defined(PETSC_HAVE_PASTIX) 4146 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 4147 #endif 4148 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 4149 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 4150 #endif 4151 #if defined(PETSC_HAVE_SUPERLU) 4152 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 4153 #endif 4154 #if defined(PETSC_HAVE_SUPERLU_DIST) 4155 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 4156 #endif 4157 #if defined(PETSC_HAVE_MUMPS) 4158 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_aij_mumps);CHKERRQ(ierr); 4159 #endif 4160 #if defined(PETSC_HAVE_UMFPACK) 4161 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 4162 #endif 4163 #if defined(PETSC_HAVE_CHOLMOD) 4164 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 4165 #endif 4166 #if defined(PETSC_HAVE_LUSOL) 4167 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 4168 #endif 4169 #if defined(PETSC_HAVE_CLIQUE) 4170 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C",MatGetFactor_aij_clique);CHKERRQ(ierr); 4171 #endif 4172 4173 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 4174 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 4175 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 4176 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4177 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4178 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4179 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4180 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4181 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4182 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4183 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4184 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4185 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4186 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4187 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4188 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4189 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4190 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4191 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4192 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4193 PetscFunctionReturn(0); 4194 } 4195 4196 #undef __FUNCT__ 4197 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4198 /* 4199 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4200 */ 4201 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4202 { 4203 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4204 PetscErrorCode ierr; 4205 PetscInt i,m = A->rmap->n; 4206 4207 PetscFunctionBegin; 4208 c = (Mat_SeqAIJ*)C->data; 4209 4210 C->factortype = A->factortype; 4211 c->row = 0; 4212 c->col = 0; 4213 c->icol = 0; 4214 c->reallocs = 0; 4215 4216 C->assembled = PETSC_TRUE; 4217 4218 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4219 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4220 4221 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4222 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4223 for (i=0; i<m; i++) { 4224 c->imax[i] = a->imax[i]; 4225 c->ilen[i] = a->ilen[i]; 4226 } 4227 4228 /* allocate the matrix space */ 4229 if (mallocmatspace) { 4230 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4231 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4232 4233 c->singlemalloc = PETSC_TRUE; 4234 4235 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4236 if (m > 0) { 4237 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4238 if (cpvalues == MAT_COPY_VALUES) { 4239 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4240 } else { 4241 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4242 } 4243 } 4244 } 4245 4246 c->ignorezeroentries = a->ignorezeroentries; 4247 c->roworiented = a->roworiented; 4248 c->nonew = a->nonew; 4249 if (a->diag) { 4250 ierr = PetscMalloc1((m+1),&c->diag);CHKERRQ(ierr); 4251 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4252 for (i=0; i<m; i++) { 4253 c->diag[i] = a->diag[i]; 4254 } 4255 } else c->diag = 0; 4256 4257 c->solve_work = 0; 4258 c->saved_values = 0; 4259 c->idiag = 0; 4260 c->ssor_work = 0; 4261 c->keepnonzeropattern = a->keepnonzeropattern; 4262 c->free_a = PETSC_TRUE; 4263 c->free_ij = PETSC_TRUE; 4264 c->xtoy = 0; 4265 c->XtoY = 0; 4266 4267 c->rmax = a->rmax; 4268 c->nz = a->nz; 4269 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4270 C->preallocated = PETSC_TRUE; 4271 4272 c->compressedrow.use = a->compressedrow.use; 4273 c->compressedrow.nrows = a->compressedrow.nrows; 4274 if (a->compressedrow.use) { 4275 i = a->compressedrow.nrows; 4276 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4277 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4278 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4279 } else { 4280 c->compressedrow.use = PETSC_FALSE; 4281 c->compressedrow.i = NULL; 4282 c->compressedrow.rindex = NULL; 4283 } 4284 C->same_nonzero = A->same_nonzero; 4285 4286 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4287 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4288 PetscFunctionReturn(0); 4289 } 4290 4291 #undef __FUNCT__ 4292 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4293 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4294 { 4295 PetscErrorCode ierr; 4296 4297 PetscFunctionBegin; 4298 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4299 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4300 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4301 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 4302 } 4303 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4304 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4305 PetscFunctionReturn(0); 4306 } 4307 4308 #undef __FUNCT__ 4309 #define __FUNCT__ "MatLoad_SeqAIJ" 4310 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4311 { 4312 Mat_SeqAIJ *a; 4313 PetscErrorCode ierr; 4314 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4315 int fd; 4316 PetscMPIInt size; 4317 MPI_Comm comm; 4318 PetscInt bs = 1; 4319 4320 PetscFunctionBegin; 4321 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4322 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4323 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4324 4325 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4326 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4327 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4328 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4329 4330 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4331 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4332 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4333 M = header[1]; N = header[2]; nz = header[3]; 4334 4335 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4336 4337 /* read in row lengths */ 4338 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4339 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4340 4341 /* check if sum of rowlengths is same as nz */ 4342 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4343 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 4344 4345 /* set global size if not set already*/ 4346 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4347 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4348 } else { 4349 /* if sizes and type are already set, check if the vector global sizes are correct */ 4350 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4351 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4352 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4353 } 4354 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 4355 } 4356 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4357 a = (Mat_SeqAIJ*)newMat->data; 4358 4359 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4360 4361 /* read in nonzero values */ 4362 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4363 4364 /* set matrix "i" values */ 4365 a->i[0] = 0; 4366 for (i=1; i<= M; i++) { 4367 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4368 a->ilen[i-1] = rowlengths[i-1]; 4369 } 4370 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4371 4372 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4373 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4374 PetscFunctionReturn(0); 4375 } 4376 4377 #undef __FUNCT__ 4378 #define __FUNCT__ "MatEqual_SeqAIJ" 4379 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4380 { 4381 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4382 PetscErrorCode ierr; 4383 #if defined(PETSC_USE_COMPLEX) 4384 PetscInt k; 4385 #endif 4386 4387 PetscFunctionBegin; 4388 /* If the matrix dimensions are not equal,or no of nonzeros */ 4389 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4390 *flg = PETSC_FALSE; 4391 PetscFunctionReturn(0); 4392 } 4393 4394 /* if the a->i are the same */ 4395 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4396 if (!*flg) PetscFunctionReturn(0); 4397 4398 /* if a->j are the same */ 4399 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4400 if (!*flg) PetscFunctionReturn(0); 4401 4402 /* if a->a are the same */ 4403 #if defined(PETSC_USE_COMPLEX) 4404 for (k=0; k<a->nz; k++) { 4405 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4406 *flg = PETSC_FALSE; 4407 PetscFunctionReturn(0); 4408 } 4409 } 4410 #else 4411 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4412 #endif 4413 PetscFunctionReturn(0); 4414 } 4415 4416 #undef __FUNCT__ 4417 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4418 /*@ 4419 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4420 provided by the user. 4421 4422 Collective on MPI_Comm 4423 4424 Input Parameters: 4425 + comm - must be an MPI communicator of size 1 4426 . m - number of rows 4427 . n - number of columns 4428 . i - row indices 4429 . j - column indices 4430 - a - matrix values 4431 4432 Output Parameter: 4433 . mat - the matrix 4434 4435 Level: intermediate 4436 4437 Notes: 4438 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4439 once the matrix is destroyed and not before 4440 4441 You cannot set new nonzero locations into this matrix, that will generate an error. 4442 4443 The i and j indices are 0 based 4444 4445 The format which is used for the sparse matrix input, is equivalent to a 4446 row-major ordering.. i.e for the following matrix, the input data expected is 4447 as shown: 4448 4449 1 0 0 4450 2 0 3 4451 4 5 6 4452 4453 i = {0,1,3,6} [size = nrow+1 = 3+1] 4454 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4455 v = {1,2,3,4,5,6} [size = nz = 6] 4456 4457 4458 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4459 4460 @*/ 4461 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4462 { 4463 PetscErrorCode ierr; 4464 PetscInt ii; 4465 Mat_SeqAIJ *aij; 4466 #if defined(PETSC_USE_DEBUG) 4467 PetscInt jj; 4468 #endif 4469 4470 PetscFunctionBegin; 4471 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4472 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4473 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4474 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4475 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4476 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4477 aij = (Mat_SeqAIJ*)(*mat)->data; 4478 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4479 4480 aij->i = i; 4481 aij->j = j; 4482 aij->a = a; 4483 aij->singlemalloc = PETSC_FALSE; 4484 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4485 aij->free_a = PETSC_FALSE; 4486 aij->free_ij = PETSC_FALSE; 4487 4488 for (ii=0; ii<m; ii++) { 4489 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4490 #if defined(PETSC_USE_DEBUG) 4491 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 4492 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4493 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 4494 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 4495 } 4496 #endif 4497 } 4498 #if defined(PETSC_USE_DEBUG) 4499 for (ii=0; ii<aij->i[m]; ii++) { 4500 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4501 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 4502 } 4503 #endif 4504 4505 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4506 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4507 PetscFunctionReturn(0); 4508 } 4509 #undef __FUNCT__ 4510 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4511 /*@C 4512 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4513 provided by the user. 4514 4515 Collective on MPI_Comm 4516 4517 Input Parameters: 4518 + comm - must be an MPI communicator of size 1 4519 . m - number of rows 4520 . n - number of columns 4521 . i - row indices 4522 . j - column indices 4523 . a - matrix values 4524 . nz - number of nonzeros 4525 - idx - 0 or 1 based 4526 4527 Output Parameter: 4528 . mat - the matrix 4529 4530 Level: intermediate 4531 4532 Notes: 4533 The i and j indices are 0 based 4534 4535 The format which is used for the sparse matrix input, is equivalent to a 4536 row-major ordering.. i.e for the following matrix, the input data expected is 4537 as shown: 4538 4539 1 0 0 4540 2 0 3 4541 4 5 6 4542 4543 i = {0,1,1,2,2,2} 4544 j = {0,0,2,0,1,2} 4545 v = {1,2,3,4,5,6} 4546 4547 4548 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4549 4550 @*/ 4551 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4552 { 4553 PetscErrorCode ierr; 4554 PetscInt ii, *nnz, one = 1,row,col; 4555 4556 4557 PetscFunctionBegin; 4558 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4559 for (ii = 0; ii < nz; ii++) { 4560 nnz[i[ii] - !!idx] += 1; 4561 } 4562 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4563 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4564 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4565 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4566 for (ii = 0; ii < nz; ii++) { 4567 if (idx) { 4568 row = i[ii] - 1; 4569 col = j[ii] - 1; 4570 } else { 4571 row = i[ii]; 4572 col = j[ii]; 4573 } 4574 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4575 } 4576 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4577 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4578 ierr = PetscFree(nnz);CHKERRQ(ierr); 4579 PetscFunctionReturn(0); 4580 } 4581 4582 #undef __FUNCT__ 4583 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4584 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4585 { 4586 PetscErrorCode ierr; 4587 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4588 4589 PetscFunctionBegin; 4590 if (coloring->ctype == IS_COLORING_GLOBAL) { 4591 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4592 a->coloring = coloring; 4593 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4594 PetscInt i,*larray; 4595 ISColoring ocoloring; 4596 ISColoringValue *colors; 4597 4598 /* set coloring for diagonal portion */ 4599 ierr = PetscMalloc1(A->cmap->n,&larray);CHKERRQ(ierr); 4600 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4601 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4602 ierr = PetscMalloc1(A->cmap->n,&colors);CHKERRQ(ierr); 4603 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4604 ierr = PetscFree(larray);CHKERRQ(ierr); 4605 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4606 a->coloring = ocoloring; 4607 } 4608 PetscFunctionReturn(0); 4609 } 4610 4611 #undef __FUNCT__ 4612 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4613 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4614 { 4615 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4616 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4617 MatScalar *v = a->a; 4618 PetscScalar *values = (PetscScalar*)advalues; 4619 ISColoringValue *color; 4620 4621 PetscFunctionBegin; 4622 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4623 color = a->coloring->colors; 4624 /* loop over rows */ 4625 for (i=0; i<m; i++) { 4626 nz = ii[i+1] - ii[i]; 4627 /* loop over columns putting computed value into matrix */ 4628 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4629 values += nl; /* jump to next row of derivatives */ 4630 } 4631 PetscFunctionReturn(0); 4632 } 4633 4634 #undef __FUNCT__ 4635 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4636 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4637 { 4638 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4639 PetscErrorCode ierr; 4640 4641 PetscFunctionBegin; 4642 a->idiagvalid = PETSC_FALSE; 4643 a->ibdiagvalid = PETSC_FALSE; 4644 4645 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4646 PetscFunctionReturn(0); 4647 } 4648 4649 /* 4650 Special version for direct calls from Fortran 4651 */ 4652 #include <petsc-private/fortranimpl.h> 4653 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4654 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4655 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4656 #define matsetvaluesseqaij_ matsetvaluesseqaij 4657 #endif 4658 4659 /* Change these macros so can be used in void function */ 4660 #undef CHKERRQ 4661 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4662 #undef SETERRQ2 4663 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4664 #undef SETERRQ3 4665 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4666 4667 #undef __FUNCT__ 4668 #define __FUNCT__ "matsetvaluesseqaij_" 4669 PETSC_EXTERN void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4670 { 4671 Mat A = *AA; 4672 PetscInt m = *mm, n = *nn; 4673 InsertMode is = *isis; 4674 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4675 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4676 PetscInt *imax,*ai,*ailen; 4677 PetscErrorCode ierr; 4678 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4679 MatScalar *ap,value,*aa; 4680 PetscBool ignorezeroentries = a->ignorezeroentries; 4681 PetscBool roworiented = a->roworiented; 4682 4683 PetscFunctionBegin; 4684 MatCheckPreallocated(A,1); 4685 imax = a->imax; 4686 ai = a->i; 4687 ailen = a->ilen; 4688 aj = a->j; 4689 aa = a->a; 4690 4691 for (k=0; k<m; k++) { /* loop over added rows */ 4692 row = im[k]; 4693 if (row < 0) continue; 4694 #if defined(PETSC_USE_DEBUG) 4695 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4696 #endif 4697 rp = aj + ai[row]; ap = aa + ai[row]; 4698 rmax = imax[row]; nrow = ailen[row]; 4699 low = 0; 4700 high = nrow; 4701 for (l=0; l<n; l++) { /* loop over added columns */ 4702 if (in[l] < 0) continue; 4703 #if defined(PETSC_USE_DEBUG) 4704 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4705 #endif 4706 col = in[l]; 4707 if (roworiented) value = v[l + k*n]; 4708 else value = v[k + l*m]; 4709 4710 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4711 4712 if (col <= lastcol) low = 0; 4713 else high = nrow; 4714 lastcol = col; 4715 while (high-low > 5) { 4716 t = (low+high)/2; 4717 if (rp[t] > col) high = t; 4718 else low = t; 4719 } 4720 for (i=low; i<high; i++) { 4721 if (rp[i] > col) break; 4722 if (rp[i] == col) { 4723 if (is == ADD_VALUES) ap[i] += value; 4724 else ap[i] = value; 4725 goto noinsert; 4726 } 4727 } 4728 if (value == 0.0 && ignorezeroentries) goto noinsert; 4729 if (nonew == 1) goto noinsert; 4730 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4731 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4732 N = nrow++ - 1; a->nz++; high++; 4733 /* shift up all the later entries in this row */ 4734 for (ii=N; ii>=i; ii--) { 4735 rp[ii+1] = rp[ii]; 4736 ap[ii+1] = ap[ii]; 4737 } 4738 rp[i] = col; 4739 ap[i] = value; 4740 noinsert:; 4741 low = i + 1; 4742 } 4743 ailen[row] = nrow; 4744 } 4745 A->same_nonzero = PETSC_FALSE; 4746 PetscFunctionReturnVoid(); 4747 } 4748 4749 4750