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