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