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