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