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