1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 8 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 9 #include <petscblaslapack.h> 10 #include <petscbt.h> 11 #include <petsc-private/kernels/blocktranspose.h> 12 #if defined(PETSC_THREADCOMM_ACTIVE) 13 #include <petscthreadcomm.h> 14 #endif 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 18 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 19 { 20 PetscErrorCode ierr; 21 PetscInt i,m,n; 22 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 23 24 PetscFunctionBegin; 25 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 26 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 27 if (type == NORM_2) { 28 for (i=0; i<aij->i[m]; i++) { 29 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 30 } 31 } else if (type == NORM_1) { 32 for (i=0; i<aij->i[m]; i++) { 33 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 34 } 35 } else if (type == NORM_INFINITY) { 36 for (i=0; i<aij->i[m]; i++) { 37 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 38 } 39 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 40 41 if (type == NORM_2) { 42 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ_Private" 49 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 50 { 51 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 52 const MatScalar *aa = a->a; 53 PetscInt i,m=A->rmap->n,cnt = 0; 54 const PetscInt *jj = a->j,*diag; 55 PetscInt *rows; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 60 diag = a->diag; 61 for (i=0; i<m; i++) { 62 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 63 cnt++; 64 } 65 } 66 ierr = PetscMalloc(cnt*sizeof(PetscInt),&rows);CHKERRQ(ierr); 67 cnt = 0; 68 for (i=0; i<m; i++) { 69 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 70 rows[cnt++] = i; 71 } 72 } 73 *nrows = cnt; 74 *zrows = rows; 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ" 80 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 81 { 82 PetscInt nrows,*rows; 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 *zrows = NULL; 87 ierr = MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);CHKERRQ(ierr); 88 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatFindNonzeroRows_SeqAIJ" 94 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 95 { 96 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 97 const MatScalar *aa; 98 PetscInt m=A->rmap->n,cnt = 0; 99 const PetscInt *ii; 100 PetscInt n,i,j,*rows; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 *keptrows = 0; 105 ii = a->i; 106 for (i=0; i<m; i++) { 107 n = ii[i+1] - ii[i]; 108 if (!n) { 109 cnt++; 110 goto ok1; 111 } 112 aa = a->a + ii[i]; 113 for (j=0; j<n; j++) { 114 if (aa[j] != 0.0) goto ok1; 115 } 116 cnt++; 117 ok1:; 118 } 119 if (!cnt) PetscFunctionReturn(0); 120 ierr = PetscMalloc((A->rmap->n-cnt)*sizeof(PetscInt),&rows);CHKERRQ(ierr); 121 cnt = 0; 122 for (i=0; i<m; i++) { 123 n = ii[i+1] - ii[i]; 124 if (!n) continue; 125 aa = a->a + ii[i]; 126 for (j=0; j<n; j++) { 127 if (aa[j] != 0.0) { 128 rows[cnt++] = i; 129 break; 130 } 131 } 132 } 133 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 139 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 140 { 141 PetscErrorCode ierr; 142 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 143 PetscInt i,*diag, m = Y->rmap->n; 144 MatScalar *aa = aij->a; 145 PetscScalar *v; 146 PetscBool missing; 147 148 PetscFunctionBegin; 149 if (Y->assembled) { 150 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);CHKERRQ(ierr); 151 if (!missing) { 152 diag = aij->diag; 153 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 154 if (is == INSERT_VALUES) { 155 for (i=0; i<m; i++) { 156 aa[diag[i]] = v[i]; 157 } 158 } else { 159 for (i=0; i<m; i++) { 160 aa[diag[i]] += v[i]; 161 } 162 } 163 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 167 } 168 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 174 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 175 { 176 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 177 PetscErrorCode ierr; 178 PetscInt i,ishift; 179 180 PetscFunctionBegin; 181 *m = A->rmap->n; 182 if (!ia) PetscFunctionReturn(0); 183 ishift = 0; 184 if (symmetric && !A->structurally_symmetric) { 185 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 186 } else if (oshift == 1) { 187 PetscInt *tia; 188 PetscInt nz = a->i[A->rmap->n]; 189 /* malloc space and add 1 to i and j indices */ 190 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),&tia);CHKERRQ(ierr); 191 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 192 *ia = tia; 193 if (ja) { 194 PetscInt *tja; 195 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&tja);CHKERRQ(ierr); 196 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 197 *ja = tja; 198 } 199 } else { 200 *ia = a->i; 201 if (ja) *ja = a->j; 202 } 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 208 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 209 { 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 if (!ia) PetscFunctionReturn(0); 214 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 215 ierr = PetscFree(*ia);CHKERRQ(ierr); 216 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 217 } 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 223 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 224 { 225 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 226 PetscErrorCode ierr; 227 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 228 PetscInt nz = a->i[m],row,*jj,mr,col; 229 230 PetscFunctionBegin; 231 *nn = n; 232 if (!ia) PetscFunctionReturn(0); 233 if (symmetric) { 234 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 235 } else { 236 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 237 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 238 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 239 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 240 jj = a->j; 241 for (i=0; i<nz; i++) { 242 collengths[jj[i]]++; 243 } 244 cia[0] = oshift; 245 for (i=0; i<n; i++) { 246 cia[i+1] = cia[i] + collengths[i]; 247 } 248 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 249 jj = a->j; 250 for (row=0; row<m; row++) { 251 mr = a->i[row+1] - a->i[row]; 252 for (i=0; i<mr; i++) { 253 col = *jj++; 254 255 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 256 } 257 } 258 ierr = PetscFree(collengths);CHKERRQ(ierr); 259 *ia = cia; *ja = cja; 260 } 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 266 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 267 { 268 PetscErrorCode ierr; 269 270 PetscFunctionBegin; 271 if (!ia) PetscFunctionReturn(0); 272 273 ierr = PetscFree(*ia);CHKERRQ(ierr); 274 ierr = PetscFree(*ja);CHKERRQ(ierr); 275 PetscFunctionReturn(0); 276 } 277 278 #undef __FUNCT__ 279 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 280 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 281 { 282 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 283 PetscInt *ai = a->i; 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatSetValues_SeqAIJ" 293 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 294 { 295 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 296 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 297 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 298 PetscErrorCode ierr; 299 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 300 MatScalar *ap,value,*aa = a->a; 301 PetscBool ignorezeroentries = a->ignorezeroentries; 302 PetscBool roworiented = a->roworiented; 303 304 PetscFunctionBegin; 305 if (v) PetscValidScalarPointer(v,6); 306 for (k=0; k<m; k++) { /* loop over added rows */ 307 row = im[k]; 308 if (row < 0) continue; 309 #if defined(PETSC_USE_DEBUG) 310 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 311 #endif 312 rp = aj + ai[row]; ap = aa + ai[row]; 313 rmax = imax[row]; nrow = ailen[row]; 314 low = 0; 315 high = nrow; 316 for (l=0; l<n; l++) { /* loop over added columns */ 317 if (in[l] < 0) continue; 318 #if defined(PETSC_USE_DEBUG) 319 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 320 #endif 321 col = in[l]; 322 if (v) { 323 if (roworiented) { 324 value = v[l + k*n]; 325 } else { 326 value = v[k + l*m]; 327 } 328 } else { 329 value = 0.; 330 } 331 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 332 333 if (col <= lastcol) low = 0; 334 else high = nrow; 335 lastcol = col; 336 while (high-low > 5) { 337 t = (low+high)/2; 338 if (rp[t] > col) high = t; 339 else low = t; 340 } 341 for (i=low; i<high; i++) { 342 if (rp[i] > col) break; 343 if (rp[i] == col) { 344 if (is == ADD_VALUES) ap[i] += value; 345 else ap[i] = value; 346 low = i + 1; 347 goto noinsert; 348 } 349 } 350 if (value == 0.0 && ignorezeroentries) goto noinsert; 351 if (nonew == 1) goto noinsert; 352 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 353 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 354 N = nrow++ - 1; a->nz++; high++; 355 /* shift up all the later entries in this row */ 356 for (ii=N; ii>=i; ii--) { 357 rp[ii+1] = rp[ii]; 358 ap[ii+1] = ap[ii]; 359 } 360 rp[i] = col; 361 ap[i] = value; 362 low = i + 1; 363 noinsert:; 364 } 365 ailen[row] = nrow; 366 } 367 A->same_nonzero = PETSC_FALSE; 368 PetscFunctionReturn(0); 369 } 370 371 372 #undef __FUNCT__ 373 #define __FUNCT__ "MatGetValues_SeqAIJ" 374 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 375 { 376 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 377 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 378 PetscInt *ai = a->i,*ailen = a->ilen; 379 MatScalar *ap,*aa = a->a; 380 381 PetscFunctionBegin; 382 for (k=0; k<m; k++) { /* loop over rows */ 383 row = im[k]; 384 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 385 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 386 rp = aj + ai[row]; ap = aa + ai[row]; 387 nrow = ailen[row]; 388 for (l=0; l<n; l++) { /* loop over columns */ 389 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 390 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 391 col = in[l]; 392 high = nrow; low = 0; /* assume unsorted */ 393 while (high-low > 5) { 394 t = (low+high)/2; 395 if (rp[t] > col) high = t; 396 else low = t; 397 } 398 for (i=low; i<high; i++) { 399 if (rp[i] > col) break; 400 if (rp[i] == col) { 401 *v++ = ap[i]; 402 goto finished; 403 } 404 } 405 *v++ = 0.0; 406 finished:; 407 } 408 } 409 PetscFunctionReturn(0); 410 } 411 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "MatView_SeqAIJ_Binary" 415 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 416 { 417 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 418 PetscErrorCode ierr; 419 PetscInt i,*col_lens; 420 int fd; 421 FILE *file; 422 423 PetscFunctionBegin; 424 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 425 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 426 427 col_lens[0] = MAT_FILE_CLASSID; 428 col_lens[1] = A->rmap->n; 429 col_lens[2] = A->cmap->n; 430 col_lens[3] = a->nz; 431 432 /* store lengths of each row and write (including header) to file */ 433 for (i=0; i<A->rmap->n; i++) { 434 col_lens[4+i] = a->i[i+1] - a->i[i]; 435 } 436 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 437 ierr = PetscFree(col_lens);CHKERRQ(ierr); 438 439 /* store column indices (zero start index) */ 440 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 441 442 /* store nonzero values */ 443 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 444 445 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 446 if (file) { 447 fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs); 448 } 449 PetscFunctionReturn(0); 450 } 451 452 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 453 454 #undef __FUNCT__ 455 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 456 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 457 { 458 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 459 PetscErrorCode ierr; 460 PetscInt i,j,m = A->rmap->n,shift=0; 461 const char *name; 462 PetscViewerFormat format; 463 464 PetscFunctionBegin; 465 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 466 if (format == PETSC_VIEWER_ASCII_MATLAB) { 467 PetscInt nofinalvalue = 0; 468 if ((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 <petsc-private/kernels/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 /*139*/0, 3180 0 3181 }; 3182 3183 #undef __FUNCT__ 3184 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3185 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3186 { 3187 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3188 PetscInt i,nz,n; 3189 3190 PetscFunctionBegin; 3191 nz = aij->maxnz; 3192 n = mat->rmap->n; 3193 for (i=0; i<nz; i++) { 3194 aij->j[i] = indices[i]; 3195 } 3196 aij->nz = nz; 3197 for (i=0; i<n; i++) { 3198 aij->ilen[i] = aij->imax[i]; 3199 } 3200 PetscFunctionReturn(0); 3201 } 3202 3203 #undef __FUNCT__ 3204 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3205 /*@ 3206 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3207 in the matrix. 3208 3209 Input Parameters: 3210 + mat - the SeqAIJ matrix 3211 - indices - the column indices 3212 3213 Level: advanced 3214 3215 Notes: 3216 This can be called if you have precomputed the nonzero structure of the 3217 matrix and want to provide it to the matrix object to improve the performance 3218 of the MatSetValues() operation. 3219 3220 You MUST have set the correct numbers of nonzeros per row in the call to 3221 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3222 3223 MUST be called before any calls to MatSetValues(); 3224 3225 The indices should start with zero, not one. 3226 3227 @*/ 3228 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3229 { 3230 PetscErrorCode ierr; 3231 3232 PetscFunctionBegin; 3233 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3234 PetscValidPointer(indices,2); 3235 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3236 PetscFunctionReturn(0); 3237 } 3238 3239 /* ----------------------------------------------------------------------------------------*/ 3240 3241 #undef __FUNCT__ 3242 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3243 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3244 { 3245 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3246 PetscErrorCode ierr; 3247 size_t nz = aij->i[mat->rmap->n]; 3248 3249 PetscFunctionBegin; 3250 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3251 3252 /* allocate space for values if not already there */ 3253 if (!aij->saved_values) { 3254 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3255 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3256 } 3257 3258 /* copy values over */ 3259 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3260 PetscFunctionReturn(0); 3261 } 3262 3263 #undef __FUNCT__ 3264 #define __FUNCT__ "MatStoreValues" 3265 /*@ 3266 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3267 example, reuse of the linear part of a Jacobian, while recomputing the 3268 nonlinear portion. 3269 3270 Collect on Mat 3271 3272 Input Parameters: 3273 . mat - the matrix (currently only AIJ matrices support this option) 3274 3275 Level: advanced 3276 3277 Common Usage, with SNESSolve(): 3278 $ Create Jacobian matrix 3279 $ Set linear terms into matrix 3280 $ Apply boundary conditions to matrix, at this time matrix must have 3281 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3282 $ boundary conditions again will not change the nonzero structure 3283 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3284 $ ierr = MatStoreValues(mat); 3285 $ Call SNESSetJacobian() with matrix 3286 $ In your Jacobian routine 3287 $ ierr = MatRetrieveValues(mat); 3288 $ Set nonlinear terms in matrix 3289 3290 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3291 $ // build linear portion of Jacobian 3292 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3293 $ ierr = MatStoreValues(mat); 3294 $ loop over nonlinear iterations 3295 $ ierr = MatRetrieveValues(mat); 3296 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3297 $ // call MatAssemblyBegin/End() on matrix 3298 $ Solve linear system with Jacobian 3299 $ endloop 3300 3301 Notes: 3302 Matrix must already be assemblied before calling this routine 3303 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3304 calling this routine. 3305 3306 When this is called multiple times it overwrites the previous set of stored values 3307 and does not allocated additional space. 3308 3309 .seealso: MatRetrieveValues() 3310 3311 @*/ 3312 PetscErrorCode MatStoreValues(Mat mat) 3313 { 3314 PetscErrorCode ierr; 3315 3316 PetscFunctionBegin; 3317 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3318 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3319 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3320 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3321 PetscFunctionReturn(0); 3322 } 3323 3324 #undef __FUNCT__ 3325 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3326 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3327 { 3328 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3329 PetscErrorCode ierr; 3330 PetscInt nz = aij->i[mat->rmap->n]; 3331 3332 PetscFunctionBegin; 3333 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3334 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3335 /* copy values over */ 3336 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3337 PetscFunctionReturn(0); 3338 } 3339 3340 #undef __FUNCT__ 3341 #define __FUNCT__ "MatRetrieveValues" 3342 /*@ 3343 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3344 example, reuse of the linear part of a Jacobian, while recomputing the 3345 nonlinear portion. 3346 3347 Collect on Mat 3348 3349 Input Parameters: 3350 . mat - the matrix (currently on AIJ matrices support this option) 3351 3352 Level: advanced 3353 3354 .seealso: MatStoreValues() 3355 3356 @*/ 3357 PetscErrorCode MatRetrieveValues(Mat mat) 3358 { 3359 PetscErrorCode ierr; 3360 3361 PetscFunctionBegin; 3362 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3363 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3364 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3365 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3366 PetscFunctionReturn(0); 3367 } 3368 3369 3370 /* --------------------------------------------------------------------------------*/ 3371 #undef __FUNCT__ 3372 #define __FUNCT__ "MatCreateSeqAIJ" 3373 /*@C 3374 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3375 (the default parallel PETSc format). For good matrix assembly performance 3376 the user should preallocate the matrix storage by setting the parameter nz 3377 (or the array nnz). By setting these parameters accurately, performance 3378 during matrix assembly can be increased by more than a factor of 50. 3379 3380 Collective on MPI_Comm 3381 3382 Input Parameters: 3383 + comm - MPI communicator, set to PETSC_COMM_SELF 3384 . m - number of rows 3385 . n - number of columns 3386 . nz - number of nonzeros per row (same for all rows) 3387 - nnz - array containing the number of nonzeros in the various rows 3388 (possibly different for each row) or NULL 3389 3390 Output Parameter: 3391 . A - the matrix 3392 3393 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3394 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3395 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3396 3397 Notes: 3398 If nnz is given then nz is ignored 3399 3400 The AIJ format (also called the Yale sparse matrix format or 3401 compressed row storage), is fully compatible with standard Fortran 77 3402 storage. That is, the stored row and column indices can begin at 3403 either one (as in Fortran) or zero. See the users' manual for details. 3404 3405 Specify the preallocated storage with either nz or nnz (not both). 3406 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3407 allocation. For large problems you MUST preallocate memory or you 3408 will get TERRIBLE performance, see the users' manual chapter on matrices. 3409 3410 By default, this format uses inodes (identical nodes) when possible, to 3411 improve numerical efficiency of matrix-vector products and solves. We 3412 search for consecutive rows with the same nonzero structure, thereby 3413 reusing matrix information to achieve increased efficiency. 3414 3415 Options Database Keys: 3416 + -mat_no_inode - Do not use inodes 3417 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3418 3419 Level: intermediate 3420 3421 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3422 3423 @*/ 3424 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3425 { 3426 PetscErrorCode ierr; 3427 3428 PetscFunctionBegin; 3429 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3430 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3431 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3432 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3433 PetscFunctionReturn(0); 3434 } 3435 3436 #undef __FUNCT__ 3437 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3438 /*@C 3439 MatSeqAIJSetPreallocation - For good matrix assembly performance 3440 the user should preallocate the matrix storage by setting the parameter nz 3441 (or the array nnz). By setting these parameters accurately, performance 3442 during matrix assembly can be increased by more than a factor of 50. 3443 3444 Collective on MPI_Comm 3445 3446 Input Parameters: 3447 + B - The matrix-free 3448 . nz - number of nonzeros per row (same for all rows) 3449 - nnz - array containing the number of nonzeros in the various rows 3450 (possibly different for each row) or NULL 3451 3452 Notes: 3453 If nnz is given then nz is ignored 3454 3455 The AIJ format (also called the Yale sparse matrix format or 3456 compressed row storage), is fully compatible with standard Fortran 77 3457 storage. That is, the stored row and column indices can begin at 3458 either one (as in Fortran) or zero. See the users' manual for details. 3459 3460 Specify the preallocated storage with either nz or nnz (not both). 3461 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3462 allocation. For large problems you MUST preallocate memory or you 3463 will get TERRIBLE performance, see the users' manual chapter on matrices. 3464 3465 You can call MatGetInfo() to get information on how effective the preallocation was; 3466 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3467 You can also run with the option -info and look for messages with the string 3468 malloc in them to see if additional memory allocation was needed. 3469 3470 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3471 entries or columns indices 3472 3473 By default, this format uses inodes (identical nodes) when possible, to 3474 improve numerical efficiency of matrix-vector products and solves. We 3475 search for consecutive rows with the same nonzero structure, thereby 3476 reusing matrix information to achieve increased efficiency. 3477 3478 Options Database Keys: 3479 + -mat_no_inode - Do not use inodes 3480 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3481 - -mat_aij_oneindex - Internally use indexing starting at 1 3482 rather than 0. Note that when calling MatSetValues(), 3483 the user still MUST index entries starting at 0! 3484 3485 Level: intermediate 3486 3487 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3488 3489 @*/ 3490 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3491 { 3492 PetscErrorCode ierr; 3493 3494 PetscFunctionBegin; 3495 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3496 PetscValidType(B,1); 3497 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3498 PetscFunctionReturn(0); 3499 } 3500 3501 #undef __FUNCT__ 3502 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3503 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3504 { 3505 Mat_SeqAIJ *b; 3506 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3507 PetscErrorCode ierr; 3508 PetscInt i; 3509 3510 PetscFunctionBegin; 3511 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3512 if (nz == MAT_SKIP_ALLOCATION) { 3513 skipallocation = PETSC_TRUE; 3514 nz = 0; 3515 } 3516 3517 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3518 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3519 3520 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3521 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3522 if (nnz) { 3523 for (i=0; i<B->rmap->n; i++) { 3524 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]); 3525 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); 3526 } 3527 } 3528 3529 B->preallocated = PETSC_TRUE; 3530 3531 b = (Mat_SeqAIJ*)B->data; 3532 3533 if (!skipallocation) { 3534 if (!b->imax) { 3535 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3536 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3537 } 3538 if (!nnz) { 3539 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3540 else if (nz < 0) nz = 1; 3541 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3542 nz = nz*B->rmap->n; 3543 } else { 3544 nz = 0; 3545 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3546 } 3547 /* b->ilen will count nonzeros in each row so far. */ 3548 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3549 3550 /* allocate the matrix space */ 3551 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3552 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3553 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3554 b->i[0] = 0; 3555 for (i=1; i<B->rmap->n+1; i++) { 3556 b->i[i] = b->i[i-1] + b->imax[i-1]; 3557 } 3558 b->singlemalloc = PETSC_TRUE; 3559 b->free_a = PETSC_TRUE; 3560 b->free_ij = PETSC_TRUE; 3561 #if defined(PETSC_THREADCOMM_ACTIVE) 3562 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3563 #endif 3564 } else { 3565 b->free_a = PETSC_FALSE; 3566 b->free_ij = PETSC_FALSE; 3567 } 3568 3569 b->nz = 0; 3570 b->maxnz = nz; 3571 B->info.nz_unneeded = (double)b->maxnz; 3572 if (realalloc) { 3573 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3574 } 3575 PetscFunctionReturn(0); 3576 } 3577 3578 #undef __FUNCT__ 3579 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3580 /*@ 3581 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3582 3583 Input Parameters: 3584 + B - the matrix 3585 . i - the indices into j for the start of each row (starts with zero) 3586 . j - the column indices for each row (starts with zero) these must be sorted for each row 3587 - v - optional values in the matrix 3588 3589 Level: developer 3590 3591 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3592 3593 .keywords: matrix, aij, compressed row, sparse, sequential 3594 3595 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3596 @*/ 3597 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3598 { 3599 PetscErrorCode ierr; 3600 3601 PetscFunctionBegin; 3602 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3603 PetscValidType(B,1); 3604 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3605 PetscFunctionReturn(0); 3606 } 3607 3608 #undef __FUNCT__ 3609 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3610 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3611 { 3612 PetscInt i; 3613 PetscInt m,n; 3614 PetscInt nz; 3615 PetscInt *nnz, nz_max = 0; 3616 PetscScalar *values; 3617 PetscErrorCode ierr; 3618 3619 PetscFunctionBegin; 3620 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3621 3622 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3623 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3624 3625 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3626 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3627 for (i = 0; i < m; i++) { 3628 nz = Ii[i+1]- Ii[i]; 3629 nz_max = PetscMax(nz_max, nz); 3630 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3631 nnz[i] = nz; 3632 } 3633 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3634 ierr = PetscFree(nnz);CHKERRQ(ierr); 3635 3636 if (v) { 3637 values = (PetscScalar*) v; 3638 } else { 3639 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3640 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3641 } 3642 3643 for (i = 0; i < m; i++) { 3644 nz = Ii[i+1] - Ii[i]; 3645 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3646 } 3647 3648 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3649 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3650 3651 if (!v) { 3652 ierr = PetscFree(values);CHKERRQ(ierr); 3653 } 3654 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3655 PetscFunctionReturn(0); 3656 } 3657 3658 #include <../src/mat/impls/dense/seq/dense.h> 3659 #include <petsc-private/kernels/petscaxpy.h> 3660 3661 #undef __FUNCT__ 3662 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3663 /* 3664 Computes (B'*A')' since computing B*A directly is untenable 3665 3666 n p p 3667 ( ) ( ) ( ) 3668 m ( A ) * n ( B ) = m ( C ) 3669 ( ) ( ) ( ) 3670 3671 */ 3672 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3673 { 3674 PetscErrorCode ierr; 3675 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3676 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3677 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3678 PetscInt i,n,m,q,p; 3679 const PetscInt *ii,*idx; 3680 const PetscScalar *b,*a,*a_q; 3681 PetscScalar *c,*c_q; 3682 3683 PetscFunctionBegin; 3684 m = A->rmap->n; 3685 n = A->cmap->n; 3686 p = B->cmap->n; 3687 a = sub_a->v; 3688 b = sub_b->a; 3689 c = sub_c->v; 3690 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3691 3692 ii = sub_b->i; 3693 idx = sub_b->j; 3694 for (i=0; i<n; i++) { 3695 q = ii[i+1] - ii[i]; 3696 while (q-->0) { 3697 c_q = c + m*(*idx); 3698 a_q = a + m*i; 3699 PetscKernelAXPY(c_q,*b,a_q,m); 3700 idx++; 3701 b++; 3702 } 3703 } 3704 PetscFunctionReturn(0); 3705 } 3706 3707 #undef __FUNCT__ 3708 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3709 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3710 { 3711 PetscErrorCode ierr; 3712 PetscInt m=A->rmap->n,n=B->cmap->n; 3713 Mat Cmat; 3714 3715 PetscFunctionBegin; 3716 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); 3717 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3718 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3719 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3720 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3721 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3722 3723 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3724 3725 *C = Cmat; 3726 PetscFunctionReturn(0); 3727 } 3728 3729 /* ----------------------------------------------------------------*/ 3730 #undef __FUNCT__ 3731 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3732 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3733 { 3734 PetscErrorCode ierr; 3735 3736 PetscFunctionBegin; 3737 if (scall == MAT_INITIAL_MATRIX) { 3738 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3739 } 3740 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3741 PetscFunctionReturn(0); 3742 } 3743 3744 3745 /*MC 3746 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3747 based on compressed sparse row format. 3748 3749 Options Database Keys: 3750 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3751 3752 Level: beginner 3753 3754 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3755 M*/ 3756 3757 /*MC 3758 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3759 3760 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3761 and MATMPIAIJ otherwise. As a result, for single process communicators, 3762 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3763 for communicators controlling multiple processes. It is recommended that you call both of 3764 the above preallocation routines for simplicity. 3765 3766 Options Database Keys: 3767 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3768 3769 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3770 enough exist. 3771 3772 Level: beginner 3773 3774 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3775 M*/ 3776 3777 /*MC 3778 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3779 3780 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3781 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3782 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3783 for communicators controlling multiple processes. It is recommended that you call both of 3784 the above preallocation routines for simplicity. 3785 3786 Options Database Keys: 3787 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3788 3789 Level: beginner 3790 3791 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3792 M*/ 3793 3794 #if defined(PETSC_HAVE_PASTIX) 3795 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3796 #endif 3797 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3798 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*); 3799 #endif 3800 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3801 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3802 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3803 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*); 3804 #if defined(PETSC_HAVE_MUMPS) 3805 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3806 #endif 3807 #if defined(PETSC_HAVE_SUPERLU) 3808 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3809 #endif 3810 #if defined(PETSC_HAVE_SUPERLU_DIST) 3811 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3812 #endif 3813 #if defined(PETSC_HAVE_UMFPACK) 3814 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3815 #endif 3816 #if defined(PETSC_HAVE_CHOLMOD) 3817 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3818 #endif 3819 #if defined(PETSC_HAVE_LUSOL) 3820 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3821 #endif 3822 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3823 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3824 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3825 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3826 #endif 3827 #if defined(PETSC_HAVE_CLIQUE) 3828 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 3829 #endif 3830 3831 3832 #undef __FUNCT__ 3833 #define __FUNCT__ "MatSeqAIJGetArray" 3834 /*@C 3835 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3836 3837 Not Collective 3838 3839 Input Parameter: 3840 . mat - a MATSEQDENSE matrix 3841 3842 Output Parameter: 3843 . array - pointer to the data 3844 3845 Level: intermediate 3846 3847 .seealso: MatSeqAIJRestoreArray() 3848 @*/ 3849 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3850 { 3851 PetscErrorCode ierr; 3852 3853 PetscFunctionBegin; 3854 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3855 PetscFunctionReturn(0); 3856 } 3857 3858 #undef __FUNCT__ 3859 #define __FUNCT__ "MatSeqAIJRestoreArray" 3860 /*@C 3861 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 3862 3863 Not Collective 3864 3865 Input Parameters: 3866 . mat - a MATSEQDENSE matrix 3867 . array - pointer to the data 3868 3869 Level: intermediate 3870 3871 .seealso: MatSeqAIJGetArray() 3872 @*/ 3873 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3874 { 3875 PetscErrorCode ierr; 3876 3877 PetscFunctionBegin; 3878 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3879 PetscFunctionReturn(0); 3880 } 3881 3882 #undef __FUNCT__ 3883 #define __FUNCT__ "MatCreate_SeqAIJ" 3884 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3885 { 3886 Mat_SeqAIJ *b; 3887 PetscErrorCode ierr; 3888 PetscMPIInt size; 3889 3890 PetscFunctionBegin; 3891 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3892 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3893 3894 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3895 3896 B->data = (void*)b; 3897 3898 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3899 3900 b->row = 0; 3901 b->col = 0; 3902 b->icol = 0; 3903 b->reallocs = 0; 3904 b->ignorezeroentries = PETSC_FALSE; 3905 b->roworiented = PETSC_TRUE; 3906 b->nonew = 0; 3907 b->diag = 0; 3908 b->solve_work = 0; 3909 B->spptr = 0; 3910 b->saved_values = 0; 3911 b->idiag = 0; 3912 b->mdiag = 0; 3913 b->ssor_work = 0; 3914 b->omega = 1.0; 3915 b->fshift = 0.0; 3916 b->idiagvalid = PETSC_FALSE; 3917 b->ibdiagvalid = PETSC_FALSE; 3918 b->keepnonzeropattern = PETSC_FALSE; 3919 b->xtoy = 0; 3920 b->XtoY = 0; 3921 B->same_nonzero = PETSC_FALSE; 3922 3923 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3924 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C","MatSeqAIJGetArray_SeqAIJ",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3925 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C","MatSeqAIJRestoreArray_SeqAIJ",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3926 3927 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3928 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3929 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3930 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3931 #endif 3932 #if defined(PETSC_HAVE_PASTIX) 3933 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3934 #endif 3935 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3936 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3937 #endif 3938 #if defined(PETSC_HAVE_SUPERLU) 3939 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3940 #endif 3941 #if defined(PETSC_HAVE_SUPERLU_DIST) 3942 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3943 #endif 3944 #if defined(PETSC_HAVE_MUMPS) 3945 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3946 #endif 3947 #if defined(PETSC_HAVE_UMFPACK) 3948 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3949 #endif 3950 #if defined(PETSC_HAVE_CHOLMOD) 3951 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3952 #endif 3953 #if defined(PETSC_HAVE_LUSOL) 3954 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_aij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3955 #endif 3956 #if defined(PETSC_HAVE_CLIQUE) 3957 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C","MatGetFactor_aij_clique",MatGetFactor_aij_clique);CHKERRQ(ierr); 3958 #endif 3959 3960 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3961 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3962 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3963 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3964 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3965 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3966 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3972 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3974 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3975 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3976 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3977 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3978 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3979 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3980 PetscFunctionReturn(0); 3981 } 3982 3983 #undef __FUNCT__ 3984 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3985 /* 3986 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3987 */ 3988 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3989 { 3990 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3991 PetscErrorCode ierr; 3992 PetscInt i,m = A->rmap->n; 3993 3994 PetscFunctionBegin; 3995 c = (Mat_SeqAIJ*)C->data; 3996 3997 C->factortype = A->factortype; 3998 c->row = 0; 3999 c->col = 0; 4000 c->icol = 0; 4001 c->reallocs = 0; 4002 4003 C->assembled = PETSC_TRUE; 4004 4005 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4006 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4007 4008 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 4009 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4010 for (i=0; i<m; i++) { 4011 c->imax[i] = a->imax[i]; 4012 c->ilen[i] = a->ilen[i]; 4013 } 4014 4015 /* allocate the matrix space */ 4016 if (mallocmatspace) { 4017 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 4018 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4019 4020 c->singlemalloc = PETSC_TRUE; 4021 4022 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4023 if (m > 0) { 4024 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4025 if (cpvalues == MAT_COPY_VALUES) { 4026 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4027 } else { 4028 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4029 } 4030 } 4031 } 4032 4033 c->ignorezeroentries = a->ignorezeroentries; 4034 c->roworiented = a->roworiented; 4035 c->nonew = a->nonew; 4036 if (a->diag) { 4037 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 4038 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4039 for (i=0; i<m; i++) { 4040 c->diag[i] = a->diag[i]; 4041 } 4042 } else c->diag = 0; 4043 4044 c->solve_work = 0; 4045 c->saved_values = 0; 4046 c->idiag = 0; 4047 c->ssor_work = 0; 4048 c->keepnonzeropattern = a->keepnonzeropattern; 4049 c->free_a = PETSC_TRUE; 4050 c->free_ij = PETSC_TRUE; 4051 c->xtoy = 0; 4052 c->XtoY = 0; 4053 4054 c->rmax = a->rmax; 4055 c->nz = a->nz; 4056 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4057 C->preallocated = PETSC_TRUE; 4058 4059 c->compressedrow.use = a->compressedrow.use; 4060 c->compressedrow.nrows = a->compressedrow.nrows; 4061 c->compressedrow.check = a->compressedrow.check; 4062 if (a->compressedrow.use) { 4063 i = a->compressedrow.nrows; 4064 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 4065 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4066 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4067 } else { 4068 c->compressedrow.use = PETSC_FALSE; 4069 c->compressedrow.i = NULL; 4070 c->compressedrow.rindex = NULL; 4071 } 4072 C->same_nonzero = A->same_nonzero; 4073 4074 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4075 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4076 PetscFunctionReturn(0); 4077 } 4078 4079 #undef __FUNCT__ 4080 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4081 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4082 { 4083 PetscErrorCode ierr; 4084 4085 PetscFunctionBegin; 4086 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4087 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4088 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 4089 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4090 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4091 PetscFunctionReturn(0); 4092 } 4093 4094 #undef __FUNCT__ 4095 #define __FUNCT__ "MatLoad_SeqAIJ" 4096 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4097 { 4098 Mat_SeqAIJ *a; 4099 PetscErrorCode ierr; 4100 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4101 int fd; 4102 PetscMPIInt size; 4103 MPI_Comm comm; 4104 PetscInt bs = 1; 4105 4106 PetscFunctionBegin; 4107 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4108 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4109 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4110 4111 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4112 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4113 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4114 4115 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4116 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4117 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4118 M = header[1]; N = header[2]; nz = header[3]; 4119 4120 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4121 4122 /* read in row lengths */ 4123 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4124 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4125 4126 /* check if sum of rowlengths is same as nz */ 4127 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4128 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); 4129 4130 /* set global size if not set already*/ 4131 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4132 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4133 } else { 4134 /* if sizes and type are already set, check if the vector global sizes are correct */ 4135 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4136 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4137 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4138 } 4139 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); 4140 } 4141 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4142 a = (Mat_SeqAIJ*)newMat->data; 4143 4144 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4145 4146 /* read in nonzero values */ 4147 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4148 4149 /* set matrix "i" values */ 4150 a->i[0] = 0; 4151 for (i=1; i<= M; i++) { 4152 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4153 a->ilen[i-1] = rowlengths[i-1]; 4154 } 4155 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4156 4157 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4158 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4159 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4160 PetscFunctionReturn(0); 4161 } 4162 4163 #undef __FUNCT__ 4164 #define __FUNCT__ "MatEqual_SeqAIJ" 4165 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4166 { 4167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4168 PetscErrorCode ierr; 4169 #if defined(PETSC_USE_COMPLEX) 4170 PetscInt k; 4171 #endif 4172 4173 PetscFunctionBegin; 4174 /* If the matrix dimensions are not equal,or no of nonzeros */ 4175 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4176 *flg = PETSC_FALSE; 4177 PetscFunctionReturn(0); 4178 } 4179 4180 /* if the a->i are the same */ 4181 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4182 if (!*flg) PetscFunctionReturn(0); 4183 4184 /* if a->j are the same */ 4185 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4186 if (!*flg) PetscFunctionReturn(0); 4187 4188 /* if a->a are the same */ 4189 #if defined(PETSC_USE_COMPLEX) 4190 for (k=0; k<a->nz; k++) { 4191 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4192 *flg = PETSC_FALSE; 4193 PetscFunctionReturn(0); 4194 } 4195 } 4196 #else 4197 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4198 #endif 4199 PetscFunctionReturn(0); 4200 } 4201 4202 #undef __FUNCT__ 4203 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4204 /*@ 4205 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4206 provided by the user. 4207 4208 Collective on MPI_Comm 4209 4210 Input Parameters: 4211 + comm - must be an MPI communicator of size 1 4212 . m - number of rows 4213 . n - number of columns 4214 . i - row indices 4215 . j - column indices 4216 - a - matrix values 4217 4218 Output Parameter: 4219 . mat - the matrix 4220 4221 Level: intermediate 4222 4223 Notes: 4224 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4225 once the matrix is destroyed and not before 4226 4227 You cannot set new nonzero locations into this matrix, that will generate an error. 4228 4229 The i and j indices are 0 based 4230 4231 The format which is used for the sparse matrix input, is equivalent to a 4232 row-major ordering.. i.e for the following matrix, the input data expected is 4233 as shown: 4234 4235 1 0 0 4236 2 0 3 4237 4 5 6 4238 4239 i = {0,1,3,6} [size = nrow+1 = 3+1] 4240 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4241 v = {1,2,3,4,5,6} [size = nz = 6] 4242 4243 4244 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4245 4246 @*/ 4247 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4248 { 4249 PetscErrorCode ierr; 4250 PetscInt ii; 4251 Mat_SeqAIJ *aij; 4252 #if defined(PETSC_USE_DEBUG) 4253 PetscInt jj; 4254 #endif 4255 4256 PetscFunctionBegin; 4257 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4258 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4259 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4260 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4261 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4262 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4263 aij = (Mat_SeqAIJ*)(*mat)->data; 4264 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4265 4266 aij->i = i; 4267 aij->j = j; 4268 aij->a = a; 4269 aij->singlemalloc = PETSC_FALSE; 4270 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4271 aij->free_a = PETSC_FALSE; 4272 aij->free_ij = PETSC_FALSE; 4273 4274 for (ii=0; ii<m; ii++) { 4275 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4276 #if defined(PETSC_USE_DEBUG) 4277 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]); 4278 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4279 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); 4280 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); 4281 } 4282 #endif 4283 } 4284 #if defined(PETSC_USE_DEBUG) 4285 for (ii=0; ii<aij->i[m]; ii++) { 4286 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4287 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]); 4288 } 4289 #endif 4290 4291 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4292 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4293 PetscFunctionReturn(0); 4294 } 4295 #undef __FUNCT__ 4296 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4297 /*@C 4298 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4299 provided by the user. 4300 4301 Collective on MPI_Comm 4302 4303 Input Parameters: 4304 + comm - must be an MPI communicator of size 1 4305 . m - number of rows 4306 . n - number of columns 4307 . i - row indices 4308 . j - column indices 4309 . a - matrix values 4310 . nz - number of nonzeros 4311 - idx - 0 or 1 based 4312 4313 Output Parameter: 4314 . mat - the matrix 4315 4316 Level: intermediate 4317 4318 Notes: 4319 The i and j indices are 0 based 4320 4321 The format which is used for the sparse matrix input, is equivalent to a 4322 row-major ordering.. i.e for the following matrix, the input data expected is 4323 as shown: 4324 4325 1 0 0 4326 2 0 3 4327 4 5 6 4328 4329 i = {0,1,1,2,2,2} 4330 j = {0,0,2,0,1,2} 4331 v = {1,2,3,4,5,6} 4332 4333 4334 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4335 4336 @*/ 4337 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4338 { 4339 PetscErrorCode ierr; 4340 PetscInt ii, *nnz, one = 1,row,col; 4341 4342 4343 PetscFunctionBegin; 4344 ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 4345 ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr); 4346 for (ii = 0; ii < nz; ii++) { 4347 nnz[i[ii]] += 1; 4348 } 4349 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4350 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4351 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4352 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4353 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4354 for (ii = 0; ii < nz; ii++) { 4355 if (idx) { 4356 row = i[ii] - 1; 4357 col = j[ii] - 1; 4358 } else { 4359 row = i[ii]; 4360 col = j[ii]; 4361 } 4362 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4363 } 4364 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4365 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4366 ierr = PetscFree(nnz);CHKERRQ(ierr); 4367 PetscFunctionReturn(0); 4368 } 4369 4370 #undef __FUNCT__ 4371 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4372 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4373 { 4374 PetscErrorCode ierr; 4375 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4376 4377 PetscFunctionBegin; 4378 if (coloring->ctype == IS_COLORING_GLOBAL) { 4379 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4380 a->coloring = coloring; 4381 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4382 PetscInt i,*larray; 4383 ISColoring ocoloring; 4384 ISColoringValue *colors; 4385 4386 /* set coloring for diagonal portion */ 4387 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4388 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4389 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4390 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4391 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4392 ierr = PetscFree(larray);CHKERRQ(ierr); 4393 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4394 a->coloring = ocoloring; 4395 } 4396 PetscFunctionReturn(0); 4397 } 4398 4399 #undef __FUNCT__ 4400 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4401 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4402 { 4403 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4404 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4405 MatScalar *v = a->a; 4406 PetscScalar *values = (PetscScalar*)advalues; 4407 ISColoringValue *color; 4408 4409 PetscFunctionBegin; 4410 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4411 color = a->coloring->colors; 4412 /* loop over rows */ 4413 for (i=0; i<m; i++) { 4414 nz = ii[i+1] - ii[i]; 4415 /* loop over columns putting computed value into matrix */ 4416 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4417 values += nl; /* jump to next row of derivatives */ 4418 } 4419 PetscFunctionReturn(0); 4420 } 4421 4422 #undef __FUNCT__ 4423 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4424 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4425 { 4426 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4427 PetscErrorCode ierr; 4428 4429 PetscFunctionBegin; 4430 a->idiagvalid = PETSC_FALSE; 4431 a->ibdiagvalid = PETSC_FALSE; 4432 4433 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4434 PetscFunctionReturn(0); 4435 } 4436 4437 /* 4438 Special version for direct calls from Fortran 4439 */ 4440 #include <petsc-private/fortranimpl.h> 4441 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4442 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4443 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4444 #define matsetvaluesseqaij_ matsetvaluesseqaij 4445 #endif 4446 4447 /* Change these macros so can be used in void function */ 4448 #undef CHKERRQ 4449 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4450 #undef SETERRQ2 4451 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4452 #undef SETERRQ3 4453 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4454 4455 #undef __FUNCT__ 4456 #define __FUNCT__ "matsetvaluesseqaij_" 4457 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) 4458 { 4459 Mat A = *AA; 4460 PetscInt m = *mm, n = *nn; 4461 InsertMode is = *isis; 4462 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4463 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4464 PetscInt *imax,*ai,*ailen; 4465 PetscErrorCode ierr; 4466 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4467 MatScalar *ap,value,*aa; 4468 PetscBool ignorezeroentries = a->ignorezeroentries; 4469 PetscBool roworiented = a->roworiented; 4470 4471 PetscFunctionBegin; 4472 MatCheckPreallocated(A,1); 4473 imax = a->imax; 4474 ai = a->i; 4475 ailen = a->ilen; 4476 aj = a->j; 4477 aa = a->a; 4478 4479 for (k=0; k<m; k++) { /* loop over added rows */ 4480 row = im[k]; 4481 if (row < 0) continue; 4482 #if defined(PETSC_USE_DEBUG) 4483 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4484 #endif 4485 rp = aj + ai[row]; ap = aa + ai[row]; 4486 rmax = imax[row]; nrow = ailen[row]; 4487 low = 0; 4488 high = nrow; 4489 for (l=0; l<n; l++) { /* loop over added columns */ 4490 if (in[l] < 0) continue; 4491 #if defined(PETSC_USE_DEBUG) 4492 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4493 #endif 4494 col = in[l]; 4495 if (roworiented) value = v[l + k*n]; 4496 else value = v[k + l*m]; 4497 4498 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4499 4500 if (col <= lastcol) low = 0; 4501 else high = nrow; 4502 lastcol = col; 4503 while (high-low > 5) { 4504 t = (low+high)/2; 4505 if (rp[t] > col) high = t; 4506 else low = t; 4507 } 4508 for (i=low; i<high; i++) { 4509 if (rp[i] > col) break; 4510 if (rp[i] == col) { 4511 if (is == ADD_VALUES) ap[i] += value; 4512 else ap[i] = value; 4513 goto noinsert; 4514 } 4515 } 4516 if (value == 0.0 && ignorezeroentries) goto noinsert; 4517 if (nonew == 1) goto noinsert; 4518 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4519 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4520 N = nrow++ - 1; a->nz++; high++; 4521 /* shift up all the later entries in this row */ 4522 for (ii=N; ii>=i; ii--) { 4523 rp[ii+1] = rp[ii]; 4524 ap[ii+1] = ap[ii]; 4525 } 4526 rp[i] = col; 4527 ap[i] = value; 4528 noinsert:; 4529 low = i + 1; 4530 } 4531 ailen[row] = nrow; 4532 } 4533 A->same_nonzero = PETSC_FALSE; 4534 PetscFunctionReturnVoid(); 4535 } 4536 4537 4538