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