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