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