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