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