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