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