1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 8 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 9 #include <petscblaslapack.h> 10 #include <petscbt.h> 11 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 15 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 16 { 17 PetscErrorCode ierr; 18 PetscInt i,m,n; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 20 21 PetscFunctionBegin; 22 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 23 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 24 if (type == NORM_2) { 25 for (i=0; i<aij->i[m]; i++) { 26 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 27 } 28 } else if (type == NORM_1) { 29 for (i=0; i<aij->i[m]; i++) { 30 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 31 } 32 } else if (type == NORM_INFINITY) { 33 for (i=0; i<aij->i[m]; i++) { 34 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 35 } 36 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 37 38 if (type == NORM_2) { 39 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 40 } 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ" 46 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 47 { 48 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 49 const MatScalar *aa = a->a; 50 PetscInt i,m=A->rmap->n,cnt = 0; 51 const PetscInt *jj = a->j,*diag; 52 PetscInt *rows; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 57 diag = a->diag; 58 for (i=0; i<m; i++) { 59 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 60 cnt++; 61 } 62 } 63 ierr = PetscMalloc(cnt*sizeof(PetscInt),&rows);CHKERRQ(ierr); 64 cnt = 0; 65 for (i=0; i<m; i++) { 66 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 67 rows[cnt++] = i; 68 } 69 } 70 ierr = ISCreateGeneral(((PetscObject)A)->comm,cnt,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "MatFindNonzeroRows_SeqAIJ" 76 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 77 { 78 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 79 const MatScalar *aa; 80 PetscInt m=A->rmap->n,cnt = 0; 81 const PetscInt *ii; 82 PetscInt n,i,j,*rows; 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 *keptrows = 0; 87 ii = a->i; 88 for (i=0; i<m; i++) { 89 n = ii[i+1] - ii[i]; 90 if (!n) { 91 cnt++; 92 goto ok1; 93 } 94 aa = a->a + ii[i]; 95 for (j=0; j<n; j++) { 96 if (aa[j] != 0.0) goto ok1; 97 } 98 cnt++; 99 ok1:; 100 } 101 if (!cnt) PetscFunctionReturn(0); 102 ierr = PetscMalloc((A->rmap->n-cnt)*sizeof(PetscInt),&rows);CHKERRQ(ierr); 103 cnt = 0; 104 for (i=0; i<m; i++) { 105 n = ii[i+1] - ii[i]; 106 if (!n) continue; 107 aa = a->a + ii[i]; 108 for (j=0; j<n; j++) { 109 if (aa[j] != 0.0) { 110 rows[cnt++] = i; 111 break; 112 } 113 } 114 } 115 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 121 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 122 { 123 PetscErrorCode ierr; 124 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 125 PetscInt i,*diag, m = Y->rmap->n; 126 MatScalar *aa = aij->a; 127 PetscScalar *v; 128 PetscBool missing; 129 130 PetscFunctionBegin; 131 if (Y->assembled) { 132 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 133 if (!missing) { 134 diag = aij->diag; 135 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 136 if (is == INSERT_VALUES) { 137 for (i=0; i<m; i++) { 138 aa[diag[i]] = v[i]; 139 } 140 } else { 141 for (i=0; i<m; i++) { 142 aa[diag[i]] += v[i]; 143 } 144 } 145 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 } 149 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 155 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 156 { 157 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 158 PetscErrorCode ierr; 159 PetscInt i,ishift; 160 161 PetscFunctionBegin; 162 *m = A->rmap->n; 163 if (!ia) PetscFunctionReturn(0); 164 ishift = 0; 165 if (symmetric && !A->structurally_symmetric) { 166 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 167 } else if (oshift == 1) { 168 PetscInt nz = a->i[A->rmap->n]; 169 /* malloc space and add 1 to i and j indices */ 170 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 171 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 172 if (ja) { 173 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 174 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 175 } 176 } else { 177 *ia = a->i; 178 if (ja) *ja = a->j; 179 } 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 185 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 186 { 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 if (!ia) PetscFunctionReturn(0); 191 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 192 ierr = PetscFree(*ia);CHKERRQ(ierr); 193 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 194 } 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 200 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 201 { 202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 203 PetscErrorCode ierr; 204 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 205 PetscInt nz = a->i[m],row,*jj,mr,col; 206 207 PetscFunctionBegin; 208 *nn = n; 209 if (!ia) PetscFunctionReturn(0); 210 if (symmetric) { 211 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 212 } else { 213 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 214 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 215 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 216 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 217 jj = a->j; 218 for (i=0; i<nz; i++) { 219 collengths[jj[i]]++; 220 } 221 cia[0] = oshift; 222 for (i=0; i<n; i++) { 223 cia[i+1] = cia[i] + collengths[i]; 224 } 225 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 226 jj = a->j; 227 for (row=0; row<m; row++) { 228 mr = a->i[row+1] - a->i[row]; 229 for (i=0; i<mr; i++) { 230 col = *jj++; 231 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 232 } 233 } 234 ierr = PetscFree(collengths);CHKERRQ(ierr); 235 *ia = cia; *ja = cja; 236 } 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 242 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 243 { 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 if (!ia) PetscFunctionReturn(0); 248 249 ierr = PetscFree(*ia);CHKERRQ(ierr); 250 ierr = PetscFree(*ja);CHKERRQ(ierr); 251 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 257 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 258 { 259 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 260 PetscInt *ai = a->i; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "MatSetValues_SeqAIJ" 270 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 271 { 272 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 273 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 274 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 275 PetscErrorCode ierr; 276 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 277 MatScalar *ap,value,*aa = a->a; 278 PetscBool ignorezeroentries = a->ignorezeroentries; 279 PetscBool roworiented = a->roworiented; 280 281 PetscFunctionBegin; 282 if (v) PetscValidScalarPointer(v,6); 283 for (k=0; k<m; k++) { /* loop over added rows */ 284 row = im[k]; 285 if (row < 0) continue; 286 #if defined(PETSC_USE_DEBUG) 287 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 288 #endif 289 rp = aj + ai[row]; ap = aa + ai[row]; 290 rmax = imax[row]; nrow = ailen[row]; 291 low = 0; 292 high = nrow; 293 for (l=0; l<n; l++) { /* loop over added columns */ 294 if (in[l] < 0) continue; 295 #if defined(PETSC_USE_DEBUG) 296 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 297 #endif 298 col = in[l]; 299 if (v) { 300 if (roworiented) { 301 value = v[l + k*n]; 302 } else { 303 value = v[k + l*m]; 304 } 305 } else { 306 value = 0.; 307 } 308 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 309 310 if (col <= lastcol) low = 0; else high = nrow; 311 lastcol = col; 312 while (high-low > 5) { 313 t = (low+high)/2; 314 if (rp[t] > col) high = t; 315 else low = t; 316 } 317 for (i=low; i<high; i++) { 318 if (rp[i] > col) break; 319 if (rp[i] == col) { 320 if (is == ADD_VALUES) ap[i] += value; 321 else ap[i] = value; 322 low = i + 1; 323 goto noinsert; 324 } 325 } 326 if (value == 0.0 && ignorezeroentries) goto noinsert; 327 if (nonew == 1) goto noinsert; 328 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 329 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 330 N = nrow++ - 1; a->nz++; high++; 331 /* shift up all the later entries in this row */ 332 for (ii=N; ii>=i; ii--) { 333 rp[ii+1] = rp[ii]; 334 ap[ii+1] = ap[ii]; 335 } 336 rp[i] = col; 337 ap[i] = value; 338 low = i + 1; 339 noinsert:; 340 } 341 ailen[row] = nrow; 342 } 343 A->same_nonzero = PETSC_FALSE; 344 PetscFunctionReturn(0); 345 } 346 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "MatGetValues_SeqAIJ" 350 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 351 { 352 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 353 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 354 PetscInt *ai = a->i,*ailen = a->ilen; 355 MatScalar *ap,*aa = a->a; 356 357 PetscFunctionBegin; 358 for (k=0; k<m; k++) { /* loop over rows */ 359 row = im[k]; 360 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 361 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 362 rp = aj + ai[row]; ap = aa + ai[row]; 363 nrow = ailen[row]; 364 for (l=0; l<n; l++) { /* loop over columns */ 365 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 366 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 367 col = in[l] ; 368 high = nrow; low = 0; /* assume unsorted */ 369 while (high-low > 5) { 370 t = (low+high)/2; 371 if (rp[t] > col) high = t; 372 else low = t; 373 } 374 for (i=low; i<high; i++) { 375 if (rp[i] > col) break; 376 if (rp[i] == col) { 377 *v++ = ap[i]; 378 goto finished; 379 } 380 } 381 *v++ = 0.0; 382 finished:; 383 } 384 } 385 PetscFunctionReturn(0); 386 } 387 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "MatView_SeqAIJ_Binary" 391 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 392 { 393 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 394 PetscErrorCode ierr; 395 PetscInt i,*col_lens; 396 int fd; 397 398 PetscFunctionBegin; 399 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 400 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 401 col_lens[0] = MAT_FILE_CLASSID; 402 col_lens[1] = A->rmap->n; 403 col_lens[2] = A->cmap->n; 404 col_lens[3] = a->nz; 405 406 /* store lengths of each row and write (including header) to file */ 407 for (i=0; i<A->rmap->n; i++) { 408 col_lens[4+i] = a->i[i+1] - a->i[i]; 409 } 410 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 411 ierr = PetscFree(col_lens);CHKERRQ(ierr); 412 413 /* store column indices (zero start index) */ 414 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 415 416 /* store nonzero values */ 417 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 422 423 #undef __FUNCT__ 424 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 425 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 426 { 427 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 428 PetscErrorCode ierr; 429 PetscInt i,j,m = A->rmap->n,shift=0; 430 const char *name; 431 PetscViewerFormat format; 432 433 PetscFunctionBegin; 434 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 435 if (format == PETSC_VIEWER_ASCII_MATLAB) { 436 PetscInt nofinalvalue = 0; 437 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 438 nofinalvalue = 1; 439 } 440 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 441 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 442 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 443 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 444 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 445 446 for (i=0; i<m; i++) { 447 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 448 #if defined(PETSC_USE_COMPLEX) 449 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 450 #else 451 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 452 #endif 453 } 454 } 455 if (nofinalvalue) { 456 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 457 } 458 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 459 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 460 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 461 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 462 PetscFunctionReturn(0); 463 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 464 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 465 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 466 for (i=0; i<m; i++) { 467 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 468 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 469 #if defined(PETSC_USE_COMPLEX) 470 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 471 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 472 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 473 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 474 } else if (PetscRealPart(a->a[j]) != 0.0) { 475 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 476 } 477 #else 478 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 479 #endif 480 } 481 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 482 } 483 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 484 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 485 PetscInt nzd=0,fshift=1,*sptr; 486 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 487 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 488 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 489 for (i=0; i<m; i++) { 490 sptr[i] = nzd+1; 491 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 492 if (a->j[j] >= i) { 493 #if defined(PETSC_USE_COMPLEX) 494 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 495 #else 496 if (a->a[j] != 0.0) nzd++; 497 #endif 498 } 499 } 500 } 501 sptr[m] = nzd+1; 502 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 503 for (i=0; i<m+1; i+=6) { 504 if (i+4<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr);} 505 else if (i+3<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 506 else if (i+2<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 507 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 508 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 509 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 510 } 511 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 512 ierr = PetscFree(sptr);CHKERRQ(ierr); 513 for (i=0; i<m; i++) { 514 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 515 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 516 } 517 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 518 } 519 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 520 for (i=0; i<m; i++) { 521 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 522 if (a->j[j] >= i) { 523 #if defined(PETSC_USE_COMPLEX) 524 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 525 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 526 } 527 #else 528 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 529 #endif 530 } 531 } 532 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 533 } 534 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 535 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 536 PetscInt cnt = 0,jcnt; 537 PetscScalar value; 538 539 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 540 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 541 for (i=0; i<m; i++) { 542 jcnt = 0; 543 for (j=0; j<A->cmap->n; j++) { 544 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 545 value = a->a[cnt++]; 546 jcnt++; 547 } else { 548 value = 0.0; 549 } 550 #if defined(PETSC_USE_COMPLEX) 551 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 552 #else 553 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 554 #endif 555 } 556 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 557 } 558 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 559 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 560 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 561 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 562 #if defined(PETSC_USE_COMPLEX) 563 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 564 #else 565 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 566 #endif 567 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 568 for (i=0; i<m; i++) { 569 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 570 #if defined(PETSC_USE_COMPLEX) 571 if (PetscImaginaryPart(a->a[j]) > 0.0) { 572 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 573 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 574 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 575 } else { 576 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 577 } 578 #else 579 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 580 #endif 581 } 582 } 583 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 584 } else { 585 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 586 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 587 if (A->factortype){ 588 for (i=0; i<m; i++) { 589 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 590 /* L part */ 591 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 592 #if defined(PETSC_USE_COMPLEX) 593 if (PetscImaginaryPart(a->a[j]) > 0.0) { 594 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 595 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 596 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 597 } else { 598 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 599 } 600 #else 601 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 602 #endif 603 } 604 /* diagonal */ 605 j = a->diag[i]; 606 #if defined(PETSC_USE_COMPLEX) 607 if (PetscImaginaryPart(a->a[j]) > 0.0) { 608 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),PetscImaginaryPart(1.0/a->a[j]));CHKERRQ(ierr); 609 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 610 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),-PetscImaginaryPart(1.0/a->a[j]));CHKERRQ(ierr); 611 } else { 612 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 613 } 614 #else 615 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,1.0/a->a[j]);CHKERRQ(ierr); 616 #endif 617 618 /* U part */ 619 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 620 #if defined(PETSC_USE_COMPLEX) 621 if (PetscImaginaryPart(a->a[j]) > 0.0) { 622 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 623 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 624 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 625 } else { 626 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 627 } 628 #else 629 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 630 #endif 631 } 632 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 633 } 634 } else { 635 for (i=0; i<m; i++) { 636 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 637 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 638 #if defined(PETSC_USE_COMPLEX) 639 if (PetscImaginaryPart(a->a[j]) > 0.0) { 640 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 641 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 642 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 643 } else { 644 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 645 } 646 #else 647 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 648 #endif 649 } 650 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 651 } 652 } 653 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 654 } 655 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 659 #undef __FUNCT__ 660 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 661 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 662 { 663 Mat A = (Mat) Aa; 664 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 665 PetscErrorCode ierr; 666 PetscInt i,j,m = A->rmap->n,color; 667 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 668 PetscViewer viewer; 669 PetscViewerFormat format; 670 671 PetscFunctionBegin; 672 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 673 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 674 675 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 676 /* loop over matrix elements drawing boxes */ 677 678 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 679 /* Blue for negative, Cyan for zero and Red for positive */ 680 color = PETSC_DRAW_BLUE; 681 for (i=0; i<m; i++) { 682 y_l = m - i - 1.0; y_r = y_l + 1.0; 683 for (j=a->i[i]; j<a->i[i+1]; j++) { 684 x_l = a->j[j] ; x_r = x_l + 1.0; 685 #if defined(PETSC_USE_COMPLEX) 686 if (PetscRealPart(a->a[j]) >= 0.) continue; 687 #else 688 if (a->a[j] >= 0.) continue; 689 #endif 690 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 691 } 692 } 693 color = PETSC_DRAW_CYAN; 694 for (i=0; i<m; i++) { 695 y_l = m - i - 1.0; y_r = y_l + 1.0; 696 for (j=a->i[i]; j<a->i[i+1]; j++) { 697 x_l = a->j[j]; x_r = x_l + 1.0; 698 if (a->a[j] != 0.) continue; 699 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 700 } 701 } 702 color = PETSC_DRAW_RED; 703 for (i=0; i<m; i++) { 704 y_l = m - i - 1.0; y_r = y_l + 1.0; 705 for (j=a->i[i]; j<a->i[i+1]; j++) { 706 x_l = a->j[j]; x_r = x_l + 1.0; 707 #if defined(PETSC_USE_COMPLEX) 708 if (PetscRealPart(a->a[j]) <= 0.) continue; 709 #else 710 if (a->a[j] <= 0.) continue; 711 #endif 712 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 713 } 714 } 715 } else { 716 /* use contour shading to indicate magnitude of values */ 717 /* first determine max of all nonzero values */ 718 PetscInt nz = a->nz,count; 719 PetscDraw popup; 720 PetscReal scale; 721 722 for (i=0; i<nz; i++) { 723 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 724 } 725 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 726 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 727 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 728 count = 0; 729 for (i=0; i<m; i++) { 730 y_l = m - i - 1.0; y_r = y_l + 1.0; 731 for (j=a->i[i]; j<a->i[i+1]; j++) { 732 x_l = a->j[j]; x_r = x_l + 1.0; 733 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 734 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 735 count++; 736 } 737 } 738 } 739 PetscFunctionReturn(0); 740 } 741 742 #undef __FUNCT__ 743 #define __FUNCT__ "MatView_SeqAIJ_Draw" 744 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 745 { 746 PetscErrorCode ierr; 747 PetscDraw draw; 748 PetscReal xr,yr,xl,yl,h,w; 749 PetscBool isnull; 750 751 PetscFunctionBegin; 752 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 753 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 754 if (isnull) PetscFunctionReturn(0); 755 756 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 757 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 758 xr += w; yr += h; xl = -w; yl = -h; 759 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 760 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 761 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 762 PetscFunctionReturn(0); 763 } 764 765 #undef __FUNCT__ 766 #define __FUNCT__ "MatView_SeqAIJ" 767 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 768 { 769 PetscErrorCode ierr; 770 PetscBool iascii,isbinary,isdraw; 771 772 PetscFunctionBegin; 773 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 774 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 775 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 776 if (iascii) { 777 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 778 } else if (isbinary) { 779 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 780 } else if (isdraw) { 781 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 782 } else { 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_SeqAIJPThread" 1222 PetscErrorCode MatMult_SeqAIJPThread(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 = PetscSqrtReal(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 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 if (!a->ibdiag) { 2927 ierr = PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);CHKERRQ(ierr); 2928 ierr = PetscLogObjectMemory(A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2929 } 2930 diag = a->ibdiag; 2931 if (values) *values = a->ibdiag; 2932 /* factor and invert each block */ 2933 switch (bs){ 2934 case 1: 2935 for (i=0; i<mbs; i++) { 2936 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2937 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2938 } 2939 break; 2940 case 2: 2941 for (i=0; i<mbs; i++) { 2942 ij[0] = 2*i; ij[1] = 2*i + 1; 2943 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2944 ierr = Kernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 2945 diag += 4; 2946 } 2947 break; 2948 case 3: 2949 for (i=0; i<mbs; i++) { 2950 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2951 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2952 ierr = Kernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 2953 diag += 9; 2954 } 2955 break; 2956 case 4: 2957 for (i=0; i<mbs; i++) { 2958 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2959 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2960 ierr = Kernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 2961 diag += 16; 2962 } 2963 break; 2964 case 5: 2965 for (i=0; i<mbs; i++) { 2966 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2967 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2968 ierr = Kernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 2969 diag += 25; 2970 } 2971 break; 2972 case 6: 2973 for (i=0; i<mbs; i++) { 2974 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; 2975 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2976 ierr = Kernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 2977 diag += 36; 2978 } 2979 break; 2980 case 7: 2981 for (i=0; i<mbs; i++) { 2982 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; 2983 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 2984 ierr = Kernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 2985 diag += 49; 2986 } 2987 break; 2988 default: 2989 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 2990 for (i=0; i<mbs; i++) { 2991 for (j=0; j<bs; j++) { 2992 IJ[j] = bs*i + j; 2993 } 2994 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 2995 ierr = Kernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 2996 diag += bs2; 2997 } 2998 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 2999 } 3000 a->ibdiagvalid = PETSC_TRUE; 3001 PetscFunctionReturn(0); 3002 } 3003 3004 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 3005 /* -------------------------------------------------------------------*/ 3006 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3007 MatGetRow_SeqAIJ, 3008 MatRestoreRow_SeqAIJ, 3009 MatMult_SeqAIJ, 3010 /* 4*/ MatMultAdd_SeqAIJ, 3011 MatMultTranspose_SeqAIJ, 3012 MatMultTransposeAdd_SeqAIJ, 3013 0, 3014 0, 3015 0, 3016 /*10*/ 0, 3017 MatLUFactor_SeqAIJ, 3018 0, 3019 MatSOR_SeqAIJ, 3020 MatTranspose_SeqAIJ, 3021 /*15*/ MatGetInfo_SeqAIJ, 3022 MatEqual_SeqAIJ, 3023 MatGetDiagonal_SeqAIJ, 3024 MatDiagonalScale_SeqAIJ, 3025 MatNorm_SeqAIJ, 3026 /*20*/ 0, 3027 MatAssemblyEnd_SeqAIJ, 3028 MatSetOption_SeqAIJ, 3029 MatZeroEntries_SeqAIJ, 3030 /*24*/ MatZeroRows_SeqAIJ, 3031 0, 3032 0, 3033 0, 3034 0, 3035 /*29*/ MatSetUpPreallocation_SeqAIJ, 3036 0, 3037 0, 3038 MatGetArray_SeqAIJ, 3039 MatRestoreArray_SeqAIJ, 3040 /*34*/ MatDuplicate_SeqAIJ, 3041 0, 3042 0, 3043 MatILUFactor_SeqAIJ, 3044 0, 3045 /*39*/ MatAXPY_SeqAIJ, 3046 MatGetSubMatrices_SeqAIJ, 3047 MatIncreaseOverlap_SeqAIJ, 3048 MatGetValues_SeqAIJ, 3049 MatCopy_SeqAIJ, 3050 /*44*/ MatGetRowMax_SeqAIJ, 3051 MatScale_SeqAIJ, 3052 0, 3053 MatDiagonalSet_SeqAIJ, 3054 MatZeroRowsColumns_SeqAIJ, 3055 /*49*/ MatSetBlockSize_SeqAIJ, 3056 MatGetRowIJ_SeqAIJ, 3057 MatRestoreRowIJ_SeqAIJ, 3058 MatGetColumnIJ_SeqAIJ, 3059 MatRestoreColumnIJ_SeqAIJ, 3060 /*54*/ MatFDColoringCreate_SeqAIJ, 3061 0, 3062 0, 3063 MatPermute_SeqAIJ, 3064 0, 3065 /*59*/ 0, 3066 MatDestroy_SeqAIJ, 3067 MatView_SeqAIJ, 3068 0, 3069 0, 3070 /*64*/ 0, 3071 0, 3072 0, 3073 0, 3074 0, 3075 /*69*/ MatGetRowMaxAbs_SeqAIJ, 3076 MatGetRowMinAbs_SeqAIJ, 3077 0, 3078 MatSetColoring_SeqAIJ, 3079 #if defined(PETSC_HAVE_ADIC) 3080 MatSetValuesAdic_SeqAIJ, 3081 #else 3082 0, 3083 #endif 3084 /*74*/ MatSetValuesAdifor_SeqAIJ, 3085 MatFDColoringApply_AIJ, 3086 0, 3087 0, 3088 0, 3089 /*79*/ MatFindZeroDiagonals_SeqAIJ, 3090 0, 3091 0, 3092 0, 3093 MatLoad_SeqAIJ, 3094 /*84*/ MatIsSymmetric_SeqAIJ, 3095 MatIsHermitian_SeqAIJ, 3096 0, 3097 0, 3098 0, 3099 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 3100 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3101 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3102 MatPtAP_Basic, 3103 MatPtAPSymbolic_SeqAIJ, 3104 /*94*/ MatPtAPNumeric_SeqAIJ, 3105 MatMatMultTranspose_SeqAIJ_SeqAIJ, 3106 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 3107 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 3108 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 3109 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3110 0, 3111 0, 3112 MatConjugate_SeqAIJ, 3113 0, 3114 /*104*/MatSetValuesRow_SeqAIJ, 3115 MatRealPart_SeqAIJ, 3116 MatImaginaryPart_SeqAIJ, 3117 0, 3118 0, 3119 /*109*/MatMatSolve_SeqAIJ, 3120 0, 3121 MatGetRowMin_SeqAIJ, 3122 0, 3123 MatMissingDiagonal_SeqAIJ, 3124 /*114*/0, 3125 0, 3126 0, 3127 0, 3128 0, 3129 /*119*/0, 3130 0, 3131 0, 3132 0, 3133 MatGetMultiProcBlock_SeqAIJ, 3134 /*124*/MatFindNonzeroRows_SeqAIJ, 3135 MatGetColumnNorms_SeqAIJ, 3136 MatInvertBlockDiagonal_SeqAIJ, 3137 0, 3138 0, 3139 /*129*/0 3140 }; 3141 3142 EXTERN_C_BEGIN 3143 #undef __FUNCT__ 3144 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3145 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3146 { 3147 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3148 PetscInt i,nz,n; 3149 3150 PetscFunctionBegin; 3151 3152 nz = aij->maxnz; 3153 n = mat->rmap->n; 3154 for (i=0; i<nz; i++) { 3155 aij->j[i] = indices[i]; 3156 } 3157 aij->nz = nz; 3158 for (i=0; i<n; i++) { 3159 aij->ilen[i] = aij->imax[i]; 3160 } 3161 3162 PetscFunctionReturn(0); 3163 } 3164 EXTERN_C_END 3165 3166 #undef __FUNCT__ 3167 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3168 /*@ 3169 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3170 in the matrix. 3171 3172 Input Parameters: 3173 + mat - the SeqAIJ matrix 3174 - indices - the column indices 3175 3176 Level: advanced 3177 3178 Notes: 3179 This can be called if you have precomputed the nonzero structure of the 3180 matrix and want to provide it to the matrix object to improve the performance 3181 of the MatSetValues() operation. 3182 3183 You MUST have set the correct numbers of nonzeros per row in the call to 3184 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3185 3186 MUST be called before any calls to MatSetValues(); 3187 3188 The indices should start with zero, not one. 3189 3190 @*/ 3191 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3192 { 3193 PetscErrorCode ierr; 3194 3195 PetscFunctionBegin; 3196 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3197 PetscValidPointer(indices,2); 3198 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 3199 PetscFunctionReturn(0); 3200 } 3201 3202 /* ----------------------------------------------------------------------------------------*/ 3203 3204 EXTERN_C_BEGIN 3205 #undef __FUNCT__ 3206 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3207 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3208 { 3209 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3210 PetscErrorCode ierr; 3211 size_t nz = aij->i[mat->rmap->n]; 3212 3213 PetscFunctionBegin; 3214 if (aij->nonew != 1) { 3215 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3216 } 3217 3218 /* allocate space for values if not already there */ 3219 if (!aij->saved_values) { 3220 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3221 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3222 } 3223 3224 /* copy values over */ 3225 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3226 PetscFunctionReturn(0); 3227 } 3228 EXTERN_C_END 3229 3230 #undef __FUNCT__ 3231 #define __FUNCT__ "MatStoreValues" 3232 /*@ 3233 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3234 example, reuse of the linear part of a Jacobian, while recomputing the 3235 nonlinear portion. 3236 3237 Collect on Mat 3238 3239 Input Parameters: 3240 . mat - the matrix (currently only AIJ matrices support this option) 3241 3242 Level: advanced 3243 3244 Common Usage, with SNESSolve(): 3245 $ Create Jacobian matrix 3246 $ Set linear terms into matrix 3247 $ Apply boundary conditions to matrix, at this time matrix must have 3248 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3249 $ boundary conditions again will not change the nonzero structure 3250 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3251 $ ierr = MatStoreValues(mat); 3252 $ Call SNESSetJacobian() with matrix 3253 $ In your Jacobian routine 3254 $ ierr = MatRetrieveValues(mat); 3255 $ Set nonlinear terms in matrix 3256 3257 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3258 $ // build linear portion of Jacobian 3259 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3260 $ ierr = MatStoreValues(mat); 3261 $ loop over nonlinear iterations 3262 $ ierr = MatRetrieveValues(mat); 3263 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3264 $ // call MatAssemblyBegin/End() on matrix 3265 $ Solve linear system with Jacobian 3266 $ endloop 3267 3268 Notes: 3269 Matrix must already be assemblied before calling this routine 3270 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3271 calling this routine. 3272 3273 When this is called multiple times it overwrites the previous set of stored values 3274 and does not allocated additional space. 3275 3276 .seealso: MatRetrieveValues() 3277 3278 @*/ 3279 PetscErrorCode MatStoreValues(Mat mat) 3280 { 3281 PetscErrorCode ierr; 3282 3283 PetscFunctionBegin; 3284 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3285 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3286 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3287 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3288 PetscFunctionReturn(0); 3289 } 3290 3291 EXTERN_C_BEGIN 3292 #undef __FUNCT__ 3293 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3294 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3295 { 3296 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3297 PetscErrorCode ierr; 3298 PetscInt nz = aij->i[mat->rmap->n]; 3299 3300 PetscFunctionBegin; 3301 if (aij->nonew != 1) { 3302 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3303 } 3304 if (!aij->saved_values) { 3305 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3306 } 3307 /* copy values over */ 3308 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3309 PetscFunctionReturn(0); 3310 } 3311 EXTERN_C_END 3312 3313 #undef __FUNCT__ 3314 #define __FUNCT__ "MatRetrieveValues" 3315 /*@ 3316 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3317 example, reuse of the linear part of a Jacobian, while recomputing the 3318 nonlinear portion. 3319 3320 Collect on Mat 3321 3322 Input Parameters: 3323 . mat - the matrix (currently on AIJ matrices support this option) 3324 3325 Level: advanced 3326 3327 .seealso: MatStoreValues() 3328 3329 @*/ 3330 PetscErrorCode MatRetrieveValues(Mat mat) 3331 { 3332 PetscErrorCode ierr; 3333 3334 PetscFunctionBegin; 3335 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3336 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3337 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3338 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3339 PetscFunctionReturn(0); 3340 } 3341 3342 3343 /* --------------------------------------------------------------------------------*/ 3344 #undef __FUNCT__ 3345 #define __FUNCT__ "MatCreateSeqAIJ" 3346 /*@C 3347 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3348 (the default parallel PETSc format). For good matrix assembly performance 3349 the user should preallocate the matrix storage by setting the parameter nz 3350 (or the array nnz). By setting these parameters accurately, performance 3351 during matrix assembly can be increased by more than a factor of 50. 3352 3353 Collective on MPI_Comm 3354 3355 Input Parameters: 3356 + comm - MPI communicator, set to PETSC_COMM_SELF 3357 . m - number of rows 3358 . n - number of columns 3359 . nz - number of nonzeros per row (same for all rows) 3360 - nnz - array containing the number of nonzeros in the various rows 3361 (possibly different for each row) or PETSC_NULL 3362 3363 Output Parameter: 3364 . A - the matrix 3365 3366 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3367 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3368 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3369 3370 Notes: 3371 If nnz is given then nz is ignored 3372 3373 The AIJ format (also called the Yale sparse matrix format or 3374 compressed row storage), is fully compatible with standard Fortran 77 3375 storage. That is, the stored row and column indices can begin at 3376 either one (as in Fortran) or zero. See the users' manual for details. 3377 3378 Specify the preallocated storage with either nz or nnz (not both). 3379 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3380 allocation. For large problems you MUST preallocate memory or you 3381 will get TERRIBLE performance, see the users' manual chapter on matrices. 3382 3383 By default, this format uses inodes (identical nodes) when possible, to 3384 improve numerical efficiency of matrix-vector products and solves. We 3385 search for consecutive rows with the same nonzero structure, thereby 3386 reusing matrix information to achieve increased efficiency. 3387 3388 Options Database Keys: 3389 + -mat_no_inode - Do not use inodes 3390 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3391 3392 Level: intermediate 3393 3394 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3395 3396 @*/ 3397 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3398 { 3399 PetscErrorCode ierr; 3400 3401 PetscFunctionBegin; 3402 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3403 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3404 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3405 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3406 PetscFunctionReturn(0); 3407 } 3408 3409 #undef __FUNCT__ 3410 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3411 /*@C 3412 MatSeqAIJSetPreallocation - For good matrix assembly performance 3413 the user should preallocate the matrix storage by setting the parameter nz 3414 (or the array nnz). By setting these parameters accurately, performance 3415 during matrix assembly can be increased by more than a factor of 50. 3416 3417 Collective on MPI_Comm 3418 3419 Input Parameters: 3420 + B - The matrix-free 3421 . nz - number of nonzeros per row (same for all rows) 3422 - nnz - array containing the number of nonzeros in the various rows 3423 (possibly different for each row) or PETSC_NULL 3424 3425 Notes: 3426 If nnz is given then nz is ignored 3427 3428 The AIJ format (also called the Yale sparse matrix format or 3429 compressed row storage), is fully compatible with standard Fortran 77 3430 storage. That is, the stored row and column indices can begin at 3431 either one (as in Fortran) or zero. See the users' manual for details. 3432 3433 Specify the preallocated storage with either nz or nnz (not both). 3434 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3435 allocation. For large problems you MUST preallocate memory or you 3436 will get TERRIBLE performance, see the users' manual chapter on matrices. 3437 3438 You can call MatGetInfo() to get information on how effective the preallocation was; 3439 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3440 You can also run with the option -info and look for messages with the string 3441 malloc in them to see if additional memory allocation was needed. 3442 3443 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3444 entries or columns indices 3445 3446 By default, this format uses inodes (identical nodes) when possible, to 3447 improve numerical efficiency of matrix-vector products and solves. We 3448 search for consecutive rows with the same nonzero structure, thereby 3449 reusing matrix information to achieve increased efficiency. 3450 3451 Options Database Keys: 3452 + -mat_no_inode - Do not use inodes 3453 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3454 - -mat_aij_oneindex - Internally use indexing starting at 1 3455 rather than 0. Note that when calling MatSetValues(), 3456 the user still MUST index entries starting at 0! 3457 3458 Level: intermediate 3459 3460 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3461 3462 @*/ 3463 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3464 { 3465 PetscErrorCode ierr; 3466 3467 PetscFunctionBegin; 3468 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3469 PetscFunctionReturn(0); 3470 } 3471 3472 EXTERN_C_BEGIN 3473 #undef __FUNCT__ 3474 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3475 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3476 { 3477 Mat_SeqAIJ *b; 3478 PetscBool skipallocation = PETSC_FALSE; 3479 PetscErrorCode ierr; 3480 PetscInt i; 3481 3482 PetscFunctionBegin; 3483 3484 if (nz == MAT_SKIP_ALLOCATION) { 3485 skipallocation = PETSC_TRUE; 3486 nz = 0; 3487 } 3488 3489 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3490 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3491 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3492 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3493 3494 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3495 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3496 if (nnz) { 3497 for (i=0; i<B->rmap->n; i++) { 3498 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]); 3499 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); 3500 } 3501 } 3502 3503 B->preallocated = PETSC_TRUE; 3504 b = (Mat_SeqAIJ*)B->data; 3505 3506 if (!skipallocation) { 3507 if (!b->imax) { 3508 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3509 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3510 } 3511 if (!nnz) { 3512 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3513 else if (nz < 0) nz = 1; 3514 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3515 nz = nz*B->rmap->n; 3516 } else { 3517 nz = 0; 3518 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3519 } 3520 /* b->ilen will count nonzeros in each row so far. */ 3521 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3522 3523 /* allocate the matrix space */ 3524 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3525 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3526 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3527 b->i[0] = 0; 3528 for (i=1; i<B->rmap->n+1; i++) { 3529 b->i[i] = b->i[i-1] + b->imax[i-1]; 3530 } 3531 b->singlemalloc = PETSC_TRUE; 3532 b->free_a = PETSC_TRUE; 3533 b->free_ij = PETSC_TRUE; 3534 } else { 3535 b->free_a = PETSC_FALSE; 3536 b->free_ij = PETSC_FALSE; 3537 } 3538 3539 b->nz = 0; 3540 b->maxnz = nz; 3541 B->info.nz_unneeded = (double)b->maxnz; 3542 PetscFunctionReturn(0); 3543 } 3544 EXTERN_C_END 3545 3546 #undef __FUNCT__ 3547 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3548 /*@ 3549 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3550 3551 Input Parameters: 3552 + B - the matrix 3553 . i - the indices into j for the start of each row (starts with zero) 3554 . j - the column indices for each row (starts with zero) these must be sorted for each row 3555 - v - optional values in the matrix 3556 3557 Level: developer 3558 3559 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3560 3561 .keywords: matrix, aij, compressed row, sparse, sequential 3562 3563 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3564 @*/ 3565 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3566 { 3567 PetscErrorCode ierr; 3568 3569 PetscFunctionBegin; 3570 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3571 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3572 PetscFunctionReturn(0); 3573 } 3574 3575 EXTERN_C_BEGIN 3576 #undef __FUNCT__ 3577 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3578 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3579 { 3580 PetscInt i; 3581 PetscInt m,n; 3582 PetscInt nz; 3583 PetscInt *nnz, nz_max = 0; 3584 PetscScalar *values; 3585 PetscErrorCode ierr; 3586 3587 PetscFunctionBegin; 3588 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3589 3590 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3591 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3592 for(i = 0; i < m; i++) { 3593 nz = Ii[i+1]- Ii[i]; 3594 nz_max = PetscMax(nz_max, nz); 3595 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3596 nnz[i] = nz; 3597 } 3598 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3599 ierr = PetscFree(nnz);CHKERRQ(ierr); 3600 3601 if (v) { 3602 values = (PetscScalar*) v; 3603 } else { 3604 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3605 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3606 } 3607 3608 for(i = 0; i < m; i++) { 3609 nz = Ii[i+1] - Ii[i]; 3610 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3611 } 3612 3613 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3614 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3615 3616 if (!v) { 3617 ierr = PetscFree(values);CHKERRQ(ierr); 3618 } 3619 PetscFunctionReturn(0); 3620 } 3621 EXTERN_C_END 3622 3623 #include <../src/mat/impls/dense/seq/dense.h> 3624 #include <private/petscaxpy.h> 3625 3626 #undef __FUNCT__ 3627 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3628 /* 3629 Computes (B'*A')' since computing B*A directly is untenable 3630 3631 n p p 3632 ( ) ( ) ( ) 3633 m ( A ) * n ( B ) = m ( C ) 3634 ( ) ( ) ( ) 3635 3636 */ 3637 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3638 { 3639 PetscErrorCode ierr; 3640 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3641 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3642 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3643 PetscInt i,n,m,q,p; 3644 const PetscInt *ii,*idx; 3645 const PetscScalar *b,*a,*a_q; 3646 PetscScalar *c,*c_q; 3647 3648 PetscFunctionBegin; 3649 m = A->rmap->n; 3650 n = A->cmap->n; 3651 p = B->cmap->n; 3652 a = sub_a->v; 3653 b = sub_b->a; 3654 c = sub_c->v; 3655 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3656 3657 ii = sub_b->i; 3658 idx = sub_b->j; 3659 for (i=0; i<n; i++) { 3660 q = ii[i+1] - ii[i]; 3661 while (q-->0) { 3662 c_q = c + m*(*idx); 3663 a_q = a + m*i; 3664 PetscAXPY(c_q,*b,a_q,m); 3665 idx++; 3666 b++; 3667 } 3668 } 3669 PetscFunctionReturn(0); 3670 } 3671 3672 #undef __FUNCT__ 3673 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3674 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3675 { 3676 PetscErrorCode ierr; 3677 PetscInt m=A->rmap->n,n=B->cmap->n; 3678 Mat Cmat; 3679 3680 PetscFunctionBegin; 3681 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); 3682 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3683 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3684 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3685 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3686 Cmat->assembled = PETSC_TRUE; 3687 *C = Cmat; 3688 PetscFunctionReturn(0); 3689 } 3690 3691 /* ----------------------------------------------------------------*/ 3692 #undef __FUNCT__ 3693 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3694 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3695 { 3696 PetscErrorCode ierr; 3697 3698 PetscFunctionBegin; 3699 if (scall == MAT_INITIAL_MATRIX){ 3700 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3701 } 3702 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3703 PetscFunctionReturn(0); 3704 } 3705 3706 3707 /*MC 3708 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3709 based on compressed sparse row format. 3710 3711 Options Database Keys: 3712 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3713 3714 Level: beginner 3715 3716 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3717 M*/ 3718 3719 EXTERN_C_BEGIN 3720 #if defined(PETSC_HAVE_PASTIX) 3721 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3722 #endif 3723 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3724 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3725 #endif 3726 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3727 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3728 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3729 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3730 #if defined(PETSC_HAVE_MUMPS) 3731 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3732 #endif 3733 #if defined(PETSC_HAVE_SUPERLU) 3734 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3735 #endif 3736 #if defined(PETSC_HAVE_SUPERLU_DIST) 3737 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3738 #endif 3739 #if defined(PETSC_HAVE_SPOOLES) 3740 extern PetscErrorCode MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3741 #endif 3742 #if defined(PETSC_HAVE_UMFPACK) 3743 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3744 #endif 3745 #if defined(PETSC_HAVE_CHOLMOD) 3746 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3747 #endif 3748 #if defined(PETSC_HAVE_LUSOL) 3749 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3750 #endif 3751 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3752 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3753 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3754 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3755 #endif 3756 EXTERN_C_END 3757 3758 EXTERN_C_BEGIN 3759 #undef __FUNCT__ 3760 #define __FUNCT__ "MatCreate_SeqAIJ" 3761 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3762 { 3763 Mat_SeqAIJ *b; 3764 PetscErrorCode ierr; 3765 PetscMPIInt size; 3766 3767 PetscFunctionBegin; 3768 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3769 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3770 3771 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3772 B->data = (void*)b; 3773 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3774 b->row = 0; 3775 b->col = 0; 3776 b->icol = 0; 3777 b->reallocs = 0; 3778 b->ignorezeroentries = PETSC_FALSE; 3779 b->roworiented = PETSC_TRUE; 3780 b->nonew = 0; 3781 b->diag = 0; 3782 b->solve_work = 0; 3783 B->spptr = 0; 3784 b->saved_values = 0; 3785 b->idiag = 0; 3786 b->mdiag = 0; 3787 b->ssor_work = 0; 3788 b->omega = 1.0; 3789 b->fshift = 0.0; 3790 b->idiagvalid = PETSC_FALSE; 3791 b->ibdiagvalid = PETSC_FALSE; 3792 b->keepnonzeropattern = PETSC_FALSE; 3793 b->xtoy = 0; 3794 b->XtoY = 0; 3795 B->same_nonzero = PETSC_FALSE; 3796 3797 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3798 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3799 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3800 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3801 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3802 #endif 3803 #if defined(PETSC_HAVE_PASTIX) 3804 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3805 #endif 3806 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3807 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3808 #endif 3809 #if defined(PETSC_HAVE_SUPERLU) 3810 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3811 #endif 3812 #if defined(PETSC_HAVE_SUPERLU_DIST) 3813 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3814 #endif 3815 #if defined(PETSC_HAVE_SPOOLES) 3816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C","MatGetFactor_seqaij_spooles",MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3817 #endif 3818 #if defined(PETSC_HAVE_MUMPS) 3819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3820 #endif 3821 #if defined(PETSC_HAVE_UMFPACK) 3822 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3823 #endif 3824 #if defined(PETSC_HAVE_CHOLMOD) 3825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3826 #endif 3827 #if defined(PETSC_HAVE_LUSOL) 3828 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_seqaij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3829 #endif 3830 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3831 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3832 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3833 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3834 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3835 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3836 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3837 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3838 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3839 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3840 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3841 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3842 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3843 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3844 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3845 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3846 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3847 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3848 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3849 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3850 PetscFunctionReturn(0); 3851 } 3852 EXTERN_C_END 3853 3854 #if defined(PETSC_HAVE_PTHREADCLASSES) 3855 EXTERN_C_BEGIN 3856 #undef __FUNCT__ 3857 #define __FUNCT__ "MatCreate_SeqAIJPThread" 3858 PetscErrorCode MatCreate_SeqAIJPThread(Mat B) 3859 { 3860 PetscErrorCode ierr; 3861 3862 PetscFunctionBegin; 3863 ierr = MatCreate_SeqAIJ(B); 3864 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3865 B->ops->mult = MatMult_SeqAIJPThread; 3866 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJPTHREAD);CHKERRQ(ierr); 3867 PetscFunctionReturn(0); 3868 } 3869 EXTERN_C_END 3870 #endif 3871 3872 #undef __FUNCT__ 3873 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3874 /* 3875 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3876 */ 3877 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3878 { 3879 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3880 PetscErrorCode ierr; 3881 PetscInt i,m = A->rmap->n; 3882 3883 PetscFunctionBegin; 3884 c = (Mat_SeqAIJ*)C->data; 3885 3886 C->factortype = A->factortype; 3887 c->row = 0; 3888 c->col = 0; 3889 c->icol = 0; 3890 c->reallocs = 0; 3891 3892 C->assembled = PETSC_TRUE; 3893 3894 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3895 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3896 3897 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3898 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3899 for (i=0; i<m; i++) { 3900 c->imax[i] = a->imax[i]; 3901 c->ilen[i] = a->ilen[i]; 3902 } 3903 3904 /* allocate the matrix space */ 3905 if (mallocmatspace){ 3906 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3907 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3908 c->singlemalloc = PETSC_TRUE; 3909 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3910 if (m > 0) { 3911 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3912 if (cpvalues == MAT_COPY_VALUES) { 3913 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3914 } else { 3915 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3916 } 3917 } 3918 } 3919 3920 c->ignorezeroentries = a->ignorezeroentries; 3921 c->roworiented = a->roworiented; 3922 c->nonew = a->nonew; 3923 if (a->diag) { 3924 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3925 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3926 for (i=0; i<m; i++) { 3927 c->diag[i] = a->diag[i]; 3928 } 3929 } else c->diag = 0; 3930 c->solve_work = 0; 3931 c->saved_values = 0; 3932 c->idiag = 0; 3933 c->ssor_work = 0; 3934 c->keepnonzeropattern = a->keepnonzeropattern; 3935 c->free_a = PETSC_TRUE; 3936 c->free_ij = PETSC_TRUE; 3937 c->xtoy = 0; 3938 c->XtoY = 0; 3939 3940 c->nz = a->nz; 3941 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3942 C->preallocated = PETSC_TRUE; 3943 3944 c->compressedrow.use = a->compressedrow.use; 3945 c->compressedrow.nrows = a->compressedrow.nrows; 3946 c->compressedrow.check = a->compressedrow.check; 3947 if (a->compressedrow.use){ 3948 i = a->compressedrow.nrows; 3949 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3950 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3951 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3952 } else { 3953 c->compressedrow.use = PETSC_FALSE; 3954 c->compressedrow.i = PETSC_NULL; 3955 c->compressedrow.rindex = PETSC_NULL; 3956 } 3957 C->same_nonzero = A->same_nonzero; 3958 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3959 3960 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3961 PetscFunctionReturn(0); 3962 } 3963 3964 #undef __FUNCT__ 3965 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3966 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3967 { 3968 PetscErrorCode ierr; 3969 3970 PetscFunctionBegin; 3971 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3972 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3973 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3974 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3975 PetscFunctionReturn(0); 3976 } 3977 3978 #undef __FUNCT__ 3979 #define __FUNCT__ "MatLoad_SeqAIJ" 3980 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3981 { 3982 Mat_SeqAIJ *a; 3983 PetscErrorCode ierr; 3984 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3985 int fd; 3986 PetscMPIInt size; 3987 MPI_Comm comm; 3988 PetscInt bs = 1; 3989 3990 PetscFunctionBegin; 3991 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3992 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3993 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3994 3995 ierr = PetscOptionsBegin(comm,PETSC_NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 3996 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 3997 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3998 3999 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4000 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4001 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4002 M = header[1]; N = header[2]; nz = header[3]; 4003 4004 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4005 4006 /* read in row lengths */ 4007 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4008 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4009 4010 /* check if sum of rowlengths is same as nz */ 4011 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4012 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); 4013 4014 /* set global size if not set already*/ 4015 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4016 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4017 } else { 4018 /* if sizes and type are already set, check if the vector global sizes are correct */ 4019 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4020 if (rows < 0 && cols < 0){ /* user might provide local size instead of global size */ 4021 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4022 } 4023 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); 4024 } 4025 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4026 a = (Mat_SeqAIJ*)newMat->data; 4027 4028 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4029 4030 /* read in nonzero values */ 4031 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4032 4033 /* set matrix "i" values */ 4034 a->i[0] = 0; 4035 for (i=1; i<= M; i++) { 4036 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4037 a->ilen[i-1] = rowlengths[i-1]; 4038 } 4039 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4040 4041 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4042 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4043 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4044 PetscFunctionReturn(0); 4045 } 4046 4047 #undef __FUNCT__ 4048 #define __FUNCT__ "MatEqual_SeqAIJ" 4049 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4050 { 4051 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 4052 PetscErrorCode ierr; 4053 #if defined(PETSC_USE_COMPLEX) 4054 PetscInt k; 4055 #endif 4056 4057 PetscFunctionBegin; 4058 /* If the matrix dimensions are not equal,or no of nonzeros */ 4059 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4060 *flg = PETSC_FALSE; 4061 PetscFunctionReturn(0); 4062 } 4063 4064 /* if the a->i are the same */ 4065 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4066 if (!*flg) PetscFunctionReturn(0); 4067 4068 /* if a->j are the same */ 4069 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4070 if (!*flg) PetscFunctionReturn(0); 4071 4072 /* if a->a are the same */ 4073 #if defined(PETSC_USE_COMPLEX) 4074 for (k=0; k<a->nz; k++){ 4075 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 4076 *flg = PETSC_FALSE; 4077 PetscFunctionReturn(0); 4078 } 4079 } 4080 #else 4081 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4082 #endif 4083 PetscFunctionReturn(0); 4084 } 4085 4086 #undef __FUNCT__ 4087 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4088 /*@ 4089 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4090 provided by the user. 4091 4092 Collective on MPI_Comm 4093 4094 Input Parameters: 4095 + comm - must be an MPI communicator of size 1 4096 . m - number of rows 4097 . n - number of columns 4098 . i - row indices 4099 . j - column indices 4100 - a - matrix values 4101 4102 Output Parameter: 4103 . mat - the matrix 4104 4105 Level: intermediate 4106 4107 Notes: 4108 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4109 once the matrix is destroyed and not before 4110 4111 You cannot set new nonzero locations into this matrix, that will generate an error. 4112 4113 The i and j indices are 0 based 4114 4115 The format which is used for the sparse matrix input, is equivalent to a 4116 row-major ordering.. i.e for the following matrix, the input data expected is 4117 as shown: 4118 4119 1 0 0 4120 2 0 3 4121 4 5 6 4122 4123 i = {0,1,3,6} [size = nrow+1 = 3+1] 4124 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4125 v = {1,2,3,4,5,6} [size = nz = 6] 4126 4127 4128 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4129 4130 @*/ 4131 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 4132 { 4133 PetscErrorCode ierr; 4134 PetscInt ii; 4135 Mat_SeqAIJ *aij; 4136 #if defined(PETSC_USE_DEBUG) 4137 PetscInt jj; 4138 #endif 4139 4140 PetscFunctionBegin; 4141 if (i[0]) { 4142 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4143 } 4144 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4145 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4146 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4147 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4148 aij = (Mat_SeqAIJ*)(*mat)->data; 4149 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4150 4151 aij->i = i; 4152 aij->j = j; 4153 aij->a = a; 4154 aij->singlemalloc = PETSC_FALSE; 4155 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4156 aij->free_a = PETSC_FALSE; 4157 aij->free_ij = PETSC_FALSE; 4158 4159 for (ii=0; ii<m; ii++) { 4160 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4161 #if defined(PETSC_USE_DEBUG) 4162 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]); 4163 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4164 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); 4165 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); 4166 } 4167 #endif 4168 } 4169 #if defined(PETSC_USE_DEBUG) 4170 for (ii=0; ii<aij->i[m]; ii++) { 4171 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4172 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]); 4173 } 4174 #endif 4175 4176 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4177 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4178 PetscFunctionReturn(0); 4179 } 4180 4181 #undef __FUNCT__ 4182 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4183 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4184 { 4185 PetscErrorCode ierr; 4186 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4187 4188 PetscFunctionBegin; 4189 if (coloring->ctype == IS_COLORING_GLOBAL) { 4190 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4191 a->coloring = coloring; 4192 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4193 PetscInt i,*larray; 4194 ISColoring ocoloring; 4195 ISColoringValue *colors; 4196 4197 /* set coloring for diagonal portion */ 4198 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4199 for (i=0; i<A->cmap->n; i++) { 4200 larray[i] = i; 4201 } 4202 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 4203 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4204 for (i=0; i<A->cmap->n; i++) { 4205 colors[i] = coloring->colors[larray[i]]; 4206 } 4207 ierr = PetscFree(larray);CHKERRQ(ierr); 4208 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4209 a->coloring = ocoloring; 4210 } 4211 PetscFunctionReturn(0); 4212 } 4213 4214 #if defined(PETSC_HAVE_ADIC) 4215 EXTERN_C_BEGIN 4216 #include <adic/ad_utils.h> 4217 EXTERN_C_END 4218 4219 #undef __FUNCT__ 4220 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 4221 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 4222 { 4223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4224 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 4225 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 4226 ISColoringValue *color; 4227 4228 PetscFunctionBegin; 4229 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4230 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 4231 color = a->coloring->colors; 4232 /* loop over rows */ 4233 for (i=0; i<m; i++) { 4234 nz = ii[i+1] - ii[i]; 4235 /* loop over columns putting computed value into matrix */ 4236 for (j=0; j<nz; j++) { 4237 *v++ = values[color[*jj++]]; 4238 } 4239 values += nlen; /* jump to next row of derivatives */ 4240 } 4241 PetscFunctionReturn(0); 4242 } 4243 #endif 4244 4245 #undef __FUNCT__ 4246 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4247 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4248 { 4249 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4250 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4251 MatScalar *v = a->a; 4252 PetscScalar *values = (PetscScalar *)advalues; 4253 ISColoringValue *color; 4254 4255 PetscFunctionBegin; 4256 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4257 color = a->coloring->colors; 4258 /* loop over rows */ 4259 for (i=0; i<m; i++) { 4260 nz = ii[i+1] - ii[i]; 4261 /* loop over columns putting computed value into matrix */ 4262 for (j=0; j<nz; j++) { 4263 *v++ = values[color[*jj++]]; 4264 } 4265 values += nl; /* jump to next row of derivatives */ 4266 } 4267 PetscFunctionReturn(0); 4268 } 4269 4270 /* 4271 Special version for direct calls from Fortran 4272 */ 4273 #include <private/fortranimpl.h> 4274 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4275 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4276 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4277 #define matsetvaluesseqaij_ matsetvaluesseqaij 4278 #endif 4279 4280 /* Change these macros so can be used in void function */ 4281 #undef CHKERRQ 4282 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 4283 #undef SETERRQ2 4284 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4285 4286 EXTERN_C_BEGIN 4287 #undef __FUNCT__ 4288 #define __FUNCT__ "matsetvaluesseqaij_" 4289 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4290 { 4291 Mat A = *AA; 4292 PetscInt m = *mm, n = *nn; 4293 InsertMode is = *isis; 4294 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4295 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4296 PetscInt *imax,*ai,*ailen; 4297 PetscErrorCode ierr; 4298 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4299 MatScalar *ap,value,*aa; 4300 PetscBool ignorezeroentries = a->ignorezeroentries; 4301 PetscBool roworiented = a->roworiented; 4302 4303 PetscFunctionBegin; 4304 ierr = MatPreallocated(A);CHKERRQ(ierr); 4305 imax = a->imax; 4306 ai = a->i; 4307 ailen = a->ilen; 4308 aj = a->j; 4309 aa = a->a; 4310 4311 for (k=0; k<m; k++) { /* loop over added rows */ 4312 row = im[k]; 4313 if (row < 0) continue; 4314 #if defined(PETSC_USE_DEBUG) 4315 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4316 #endif 4317 rp = aj + ai[row]; ap = aa + ai[row]; 4318 rmax = imax[row]; nrow = ailen[row]; 4319 low = 0; 4320 high = nrow; 4321 for (l=0; l<n; l++) { /* loop over added columns */ 4322 if (in[l] < 0) continue; 4323 #if defined(PETSC_USE_DEBUG) 4324 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4325 #endif 4326 col = in[l]; 4327 if (roworiented) { 4328 value = v[l + k*n]; 4329 } else { 4330 value = v[k + l*m]; 4331 } 4332 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4333 4334 if (col <= lastcol) low = 0; else high = nrow; 4335 lastcol = col; 4336 while (high-low > 5) { 4337 t = (low+high)/2; 4338 if (rp[t] > col) high = t; 4339 else low = t; 4340 } 4341 for (i=low; i<high; i++) { 4342 if (rp[i] > col) break; 4343 if (rp[i] == col) { 4344 if (is == ADD_VALUES) ap[i] += value; 4345 else ap[i] = value; 4346 goto noinsert; 4347 } 4348 } 4349 if (value == 0.0 && ignorezeroentries) goto noinsert; 4350 if (nonew == 1) goto noinsert; 4351 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4352 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4353 N = nrow++ - 1; a->nz++; high++; 4354 /* shift up all the later entries in this row */ 4355 for (ii=N; ii>=i; ii--) { 4356 rp[ii+1] = rp[ii]; 4357 ap[ii+1] = ap[ii]; 4358 } 4359 rp[i] = col; 4360 ap[i] = value; 4361 noinsert:; 4362 low = i + 1; 4363 } 4364 ailen[row] = nrow; 4365 } 4366 A->same_nonzero = PETSC_FALSE; 4367 PetscFunctionReturnVoid(); 4368 } 4369 EXTERN_C_END 4370 4371