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