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