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