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