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