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