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