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