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