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