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