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