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_SeqPThreadAIJ" 1222 PetscErrorCode MatMult_SeqPThreadAIJ(Mat A,Vec xx,Vec yy) 1223 { 1224 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1225 PetscScalar *y; 1226 const PetscScalar *x; 1227 PetscErrorCode ierr; 1228 PetscInt m=A->rmap->n,nonzerorow=0; 1229 PetscBool usecprow=a->compressedrow.use; 1230 1231 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1232 #pragma disjoint(*x,*y,*aa) 1233 #endif 1234 1235 PetscFunctionBegin; 1236 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1237 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1238 1239 if(usecprow) { 1240 PetscInt NumPerThread,iindex; 1241 const MatScalar *aa = a->a; 1242 const PetscInt *aj = a->j,*ii = a->compressedrow.i,*ridx=a->compressedrow.rindex; 1243 PetscInt i,iStartVal,iEndVal,iStartIndex,iEndIndex; 1244 const PetscInt iNumThreads = PetscMaxThreads; //this number could be different 1245 MatMult_KernelData* kerneldatap = (MatMult_KernelData*)malloc(iNumThreads*sizeof(MatMult_KernelData)); 1246 MatMult_KernelData** pdata = (MatMult_KernelData**)malloc(iNumThreads*sizeof(MatMult_KernelData*)); 1247 1248 m = a->compressedrow.nrows; 1249 NumPerThread = ii[m]/iNumThreads; 1250 iindex = 0; 1251 for(i=0; i<iNumThreads;i++) { 1252 iStartIndex = iindex; 1253 iStartVal = ii[iStartIndex]; 1254 iEndVal = iStartVal; 1255 //determine number of rows to process 1256 while(iEndVal-iStartVal<NumPerThread) { 1257 iindex++; 1258 iEndVal = ii[iindex]; 1259 } 1260 //determine whether to go back 1 1261 if(iEndVal-iStartVal-NumPerThread>NumPerThread-(ii[iindex-1]-iStartVal)) { 1262 iindex--; 1263 iEndVal = ii[iindex]; 1264 } 1265 iEndIndex = iindex; 1266 kerneldatap[i].matdata = aa; 1267 kerneldatap[i].vecdata = x; 1268 kerneldatap[i].vecout = y; 1269 kerneldatap[i].colindnz = aj; 1270 kerneldatap[i].rownumnz = ii + iStartIndex; 1271 kerneldatap[i].numrows = iEndIndex - iStartIndex + 1; 1272 kerneldatap[i].specidx = ridx + iStartVal; 1273 kerneldatap[i].nzr = 0; 1274 pdata[i] = &kerneldatap[i]; 1275 iindex++; 1276 } 1277 ierr = MainJob(MatMult_Kernel,(void**)pdata,iNumThreads); 1278 /* collect results */ 1279 for(i=0; i<iNumThreads; i++) { 1280 nonzerorow += kerneldatap[i].nzr; 1281 } 1282 free(kerneldatap); 1283 free(pdata); 1284 } 1285 else { 1286 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1287 fortranmultaij_(&m,x,a->i,a->j,a->a,y); 1288 #else 1289 PetscInt i,iindex; 1290 const MatScalar *aa = a->a; 1291 const PetscInt *aj = a->j,*ii = a->i; 1292 const PetscInt iNumThreads = PetscMaxThreads; //this number could be different 1293 PetscInt Q = m/iNumThreads; 1294 PetscInt R = m-Q*iNumThreads; 1295 PetscBool S; 1296 1297 MatMult_KernelData* kerneldatap = (MatMult_KernelData*)malloc(iNumThreads*sizeof(MatMult_KernelData)); 1298 MatMult_KernelData** pdata = (MatMult_KernelData**)malloc(iNumThreads*sizeof(MatMult_KernelData*)); 1299 1300 iindex = 0; 1301 for(i=0; i<iNumThreads;i++) { 1302 S = (PetscBool)(i<R); 1303 kerneldatap[i].matdata = aa; 1304 kerneldatap[i].vecdata = x; 1305 kerneldatap[i].vecout = y; 1306 kerneldatap[i].colindnz = aj; 1307 kerneldatap[i].rownumnz = ii + iindex; 1308 kerneldatap[i].numrows = S?Q+1:Q; 1309 kerneldatap[i].specidx = PETSC_NULL; 1310 kerneldatap[i].nzr = iindex; //serves as the 'base' row (needed to access correctly into output vector y) 1311 pdata[i] = &kerneldatap[i]; 1312 iindex += kerneldatap[i].numrows; 1313 } 1314 MainJob(MatMult_Kernel,(void**)pdata,iNumThreads); 1315 //collect results 1316 for(i=0; i<iNumThreads; i++) { 1317 nonzerorow += kerneldatap[i].nzr; 1318 } 1319 free(kerneldatap); 1320 free(pdata); 1321 #endif 1322 } 1323 1324 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1325 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1326 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1327 PetscFunctionReturn(0); 1328 } 1329 //******************* 1330 #endif 1331 1332 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1333 #undef __FUNCT__ 1334 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1335 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1336 { 1337 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1338 PetscScalar *y,*z; 1339 const PetscScalar *x; 1340 const MatScalar *aa; 1341 PetscErrorCode ierr; 1342 PetscInt m = A->rmap->n,*aj,*ii; 1343 PetscInt n,i,*ridx=PETSC_NULL; 1344 PetscScalar sum; 1345 PetscBool usecprow=a->compressedrow.use; 1346 1347 PetscFunctionBegin; 1348 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1349 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1350 if (zz != yy) { 1351 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1352 } else { 1353 z = y; 1354 } 1355 1356 aj = a->j; 1357 aa = a->a; 1358 ii = a->i; 1359 if (usecprow){ /* use compressed row format */ 1360 if (zz != yy){ 1361 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1362 } 1363 m = a->compressedrow.nrows; 1364 ii = a->compressedrow.i; 1365 ridx = a->compressedrow.rindex; 1366 for (i=0; i<m; i++){ 1367 n = ii[i+1] - ii[i]; 1368 aj = a->j + ii[i]; 1369 aa = a->a + ii[i]; 1370 sum = y[*ridx]; 1371 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1372 z[*ridx++] = sum; 1373 } 1374 } else { /* do not use compressed row format */ 1375 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1376 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1377 #else 1378 for (i=0; i<m; i++) { 1379 n = ii[i+1] - ii[i]; 1380 aj = a->j + ii[i]; 1381 aa = a->a + ii[i]; 1382 sum = y[i]; 1383 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1384 z[i] = sum; 1385 } 1386 #endif 1387 } 1388 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1389 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1390 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1391 if (zz != yy) { 1392 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1393 } 1394 #if defined(PETSC_HAVE_CUSP) 1395 /* 1396 ierr = VecView(xx,0);CHKERRQ(ierr); 1397 ierr = VecView(zz,0);CHKERRQ(ierr); 1398 ierr = MatView(A,0);CHKERRQ(ierr); 1399 */ 1400 #endif 1401 PetscFunctionReturn(0); 1402 } 1403 1404 /* 1405 Adds diagonal pointers to sparse matrix structure. 1406 */ 1407 #undef __FUNCT__ 1408 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1409 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1410 { 1411 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1412 PetscErrorCode ierr; 1413 PetscInt i,j,m = A->rmap->n; 1414 1415 PetscFunctionBegin; 1416 if (!a->diag) { 1417 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1418 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1419 } 1420 for (i=0; i<A->rmap->n; i++) { 1421 a->diag[i] = a->i[i+1]; 1422 for (j=a->i[i]; j<a->i[i+1]; j++) { 1423 if (a->j[j] == i) { 1424 a->diag[i] = j; 1425 break; 1426 } 1427 } 1428 } 1429 PetscFunctionReturn(0); 1430 } 1431 1432 /* 1433 Checks for missing diagonals 1434 */ 1435 #undef __FUNCT__ 1436 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1437 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1438 { 1439 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1440 PetscInt *diag,*jj = a->j,i; 1441 1442 PetscFunctionBegin; 1443 *missing = PETSC_FALSE; 1444 if (A->rmap->n > 0 && !jj) { 1445 *missing = PETSC_TRUE; 1446 if (d) *d = 0; 1447 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1448 } else { 1449 diag = a->diag; 1450 for (i=0; i<A->rmap->n; i++) { 1451 if (jj[diag[i]] != i) { 1452 *missing = PETSC_TRUE; 1453 if (d) *d = i; 1454 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1455 } 1456 } 1457 } 1458 PetscFunctionReturn(0); 1459 } 1460 1461 EXTERN_C_BEGIN 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1464 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1465 { 1466 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1467 PetscErrorCode ierr; 1468 PetscInt i,*diag,m = A->rmap->n; 1469 MatScalar *v = a->a; 1470 PetscScalar *idiag,*mdiag; 1471 1472 PetscFunctionBegin; 1473 if (a->idiagvalid) PetscFunctionReturn(0); 1474 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1475 diag = a->diag; 1476 if (!a->idiag) { 1477 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1478 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1479 v = a->a; 1480 } 1481 mdiag = a->mdiag; 1482 idiag = a->idiag; 1483 1484 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1485 for (i=0; i<m; i++) { 1486 mdiag[i] = v[diag[i]]; 1487 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1488 idiag[i] = 1.0/v[diag[i]]; 1489 } 1490 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1491 } else { 1492 for (i=0; i<m; i++) { 1493 mdiag[i] = v[diag[i]]; 1494 idiag[i] = omega/(fshift + v[diag[i]]); 1495 } 1496 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1497 } 1498 a->idiagvalid = PETSC_TRUE; 1499 PetscFunctionReturn(0); 1500 } 1501 EXTERN_C_END 1502 1503 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "MatSOR_SeqAIJ" 1506 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1507 { 1508 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1509 PetscScalar *x,d,sum,*t,scale; 1510 const MatScalar *v = a->a,*idiag=0,*mdiag; 1511 const PetscScalar *b, *bs,*xb, *ts; 1512 PetscErrorCode ierr; 1513 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1514 const PetscInt *idx,*diag; 1515 1516 PetscFunctionBegin; 1517 its = its*lits; 1518 1519 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1520 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1521 a->fshift = fshift; 1522 a->omega = omega; 1523 1524 diag = a->diag; 1525 t = a->ssor_work; 1526 idiag = a->idiag; 1527 mdiag = a->mdiag; 1528 1529 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1530 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1531 CHKMEMQ; 1532 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1533 if (flag == SOR_APPLY_UPPER) { 1534 /* apply (U + D/omega) to the vector */ 1535 bs = b; 1536 for (i=0; i<m; i++) { 1537 d = fshift + mdiag[i]; 1538 n = a->i[i+1] - diag[i] - 1; 1539 idx = a->j + diag[i] + 1; 1540 v = a->a + diag[i] + 1; 1541 sum = b[i]*d/omega; 1542 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1543 x[i] = sum; 1544 } 1545 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1546 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1547 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 if (flag == SOR_APPLY_LOWER) { 1552 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1553 } else if (flag & SOR_EISENSTAT) { 1554 /* Let A = L + U + D; where L is lower trianglar, 1555 U is upper triangular, E = D/omega; This routine applies 1556 1557 (L + E)^{-1} A (U + E)^{-1} 1558 1559 to a vector efficiently using Eisenstat's trick. 1560 */ 1561 scale = (2.0/omega) - 1.0; 1562 1563 /* x = (E + U)^{-1} b */ 1564 for (i=m-1; i>=0; i--) { 1565 n = a->i[i+1] - diag[i] - 1; 1566 idx = a->j + diag[i] + 1; 1567 v = a->a + diag[i] + 1; 1568 sum = b[i]; 1569 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1570 x[i] = sum*idiag[i]; 1571 } 1572 1573 /* t = b - (2*E - D)x */ 1574 v = a->a; 1575 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1576 1577 /* t = (E + L)^{-1}t */ 1578 ts = t; 1579 diag = a->diag; 1580 for (i=0; i<m; i++) { 1581 n = diag[i] - a->i[i]; 1582 idx = a->j + a->i[i]; 1583 v = a->a + a->i[i]; 1584 sum = t[i]; 1585 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1586 t[i] = sum*idiag[i]; 1587 /* x = x + t */ 1588 x[i] += t[i]; 1589 } 1590 1591 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1592 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1593 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1594 PetscFunctionReturn(0); 1595 } 1596 if (flag & SOR_ZERO_INITIAL_GUESS) { 1597 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1598 for (i=0; i<m; i++) { 1599 n = diag[i] - a->i[i]; 1600 idx = a->j + a->i[i]; 1601 v = a->a + a->i[i]; 1602 sum = b[i]; 1603 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1604 t[i] = sum; 1605 x[i] = sum*idiag[i]; 1606 } 1607 xb = t; 1608 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1609 } else xb = b; 1610 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1611 for (i=m-1; i>=0; i--) { 1612 n = a->i[i+1] - diag[i] - 1; 1613 idx = a->j + diag[i] + 1; 1614 v = a->a + diag[i] + 1; 1615 sum = xb[i]; 1616 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1617 if (xb == b) { 1618 x[i] = sum*idiag[i]; 1619 } else { 1620 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1621 } 1622 } 1623 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1624 } 1625 its--; 1626 } 1627 while (its--) { 1628 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1629 for (i=0; i<m; i++) { 1630 n = a->i[i+1] - a->i[i]; 1631 idx = a->j + a->i[i]; 1632 v = a->a + a->i[i]; 1633 sum = b[i]; 1634 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1635 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1636 } 1637 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1638 } 1639 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1640 for (i=m-1; i>=0; i--) { 1641 n = a->i[i+1] - a->i[i]; 1642 idx = a->j + a->i[i]; 1643 v = a->a + a->i[i]; 1644 sum = b[i]; 1645 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1646 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1647 } 1648 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1649 } 1650 } 1651 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1652 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1653 CHKMEMQ; PetscFunctionReturn(0); 1654 } 1655 1656 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1659 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1660 { 1661 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1662 1663 PetscFunctionBegin; 1664 info->block_size = 1.0; 1665 info->nz_allocated = (double)a->maxnz; 1666 info->nz_used = (double)a->nz; 1667 info->nz_unneeded = (double)(a->maxnz - a->nz); 1668 info->assemblies = (double)A->num_ass; 1669 info->mallocs = (double)A->info.mallocs; 1670 info->memory = ((PetscObject)A)->mem; 1671 if (A->factortype) { 1672 info->fill_ratio_given = A->info.fill_ratio_given; 1673 info->fill_ratio_needed = A->info.fill_ratio_needed; 1674 info->factor_mallocs = A->info.factor_mallocs; 1675 } else { 1676 info->fill_ratio_given = 0; 1677 info->fill_ratio_needed = 0; 1678 info->factor_mallocs = 0; 1679 } 1680 PetscFunctionReturn(0); 1681 } 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1685 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1686 { 1687 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1688 PetscInt i,m = A->rmap->n - 1,d = 0; 1689 PetscErrorCode ierr; 1690 const PetscScalar *xx; 1691 PetscScalar *bb; 1692 PetscBool missing; 1693 1694 PetscFunctionBegin; 1695 if (x && b) { 1696 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1697 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1698 for (i=0; i<N; i++) { 1699 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1700 bb[rows[i]] = diag*xx[rows[i]]; 1701 } 1702 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1703 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1704 } 1705 1706 if (a->keepnonzeropattern) { 1707 for (i=0; i<N; i++) { 1708 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1709 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1710 } 1711 if (diag != 0.0) { 1712 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1713 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1714 for (i=0; i<N; i++) { 1715 a->a[a->diag[rows[i]]] = diag; 1716 } 1717 } 1718 A->same_nonzero = PETSC_TRUE; 1719 } else { 1720 if (diag != 0.0) { 1721 for (i=0; i<N; i++) { 1722 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1723 if (a->ilen[rows[i]] > 0) { 1724 a->ilen[rows[i]] = 1; 1725 a->a[a->i[rows[i]]] = diag; 1726 a->j[a->i[rows[i]]] = rows[i]; 1727 } else { /* in case row was completely empty */ 1728 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1729 } 1730 } 1731 } else { 1732 for (i=0; i<N; i++) { 1733 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1734 a->ilen[rows[i]] = 0; 1735 } 1736 } 1737 A->same_nonzero = PETSC_FALSE; 1738 } 1739 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1740 PetscFunctionReturn(0); 1741 } 1742 1743 #undef __FUNCT__ 1744 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1745 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1746 { 1747 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1748 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1749 PetscErrorCode ierr; 1750 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1751 const PetscScalar *xx; 1752 PetscScalar *bb; 1753 1754 PetscFunctionBegin; 1755 if (x && b) { 1756 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1757 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1758 vecs = PETSC_TRUE; 1759 } 1760 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1761 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1762 for (i=0; i<N; i++) { 1763 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1764 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1765 zeroed[rows[i]] = PETSC_TRUE; 1766 } 1767 for (i=0; i<A->rmap->n; i++) { 1768 if (!zeroed[i]) { 1769 for (j=a->i[i]; j<a->i[i+1]; j++) { 1770 if (zeroed[a->j[j]]) { 1771 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1772 a->a[j] = 0.0; 1773 } 1774 } 1775 } else if (vecs) bb[i] = diag*xx[i]; 1776 } 1777 if (x && b) { 1778 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1779 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1780 } 1781 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1782 if (diag != 0.0) { 1783 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1784 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1785 for (i=0; i<N; i++) { 1786 a->a[a->diag[rows[i]]] = diag; 1787 } 1788 } 1789 A->same_nonzero = PETSC_TRUE; 1790 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 #undef __FUNCT__ 1795 #define __FUNCT__ "MatGetRow_SeqAIJ" 1796 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1797 { 1798 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1799 PetscInt *itmp; 1800 1801 PetscFunctionBegin; 1802 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1803 1804 *nz = a->i[row+1] - a->i[row]; 1805 if (v) *v = a->a + a->i[row]; 1806 if (idx) { 1807 itmp = a->j + a->i[row]; 1808 if (*nz) { 1809 *idx = itmp; 1810 } 1811 else *idx = 0; 1812 } 1813 PetscFunctionReturn(0); 1814 } 1815 1816 /* remove this function? */ 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1819 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1820 { 1821 PetscFunctionBegin; 1822 PetscFunctionReturn(0); 1823 } 1824 1825 #undef __FUNCT__ 1826 #define __FUNCT__ "MatNorm_SeqAIJ" 1827 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1828 { 1829 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1830 MatScalar *v = a->a; 1831 PetscReal sum = 0.0; 1832 PetscErrorCode ierr; 1833 PetscInt i,j; 1834 1835 PetscFunctionBegin; 1836 if (type == NORM_FROBENIUS) { 1837 for (i=0; i<a->nz; i++) { 1838 #if defined(PETSC_USE_COMPLEX) 1839 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1840 #else 1841 sum += (*v)*(*v); v++; 1842 #endif 1843 } 1844 *nrm = 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 }; 3138 3139 EXTERN_C_BEGIN 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3142 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3143 { 3144 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3145 PetscInt i,nz,n; 3146 3147 PetscFunctionBegin; 3148 3149 nz = aij->maxnz; 3150 n = mat->rmap->n; 3151 for (i=0; i<nz; i++) { 3152 aij->j[i] = indices[i]; 3153 } 3154 aij->nz = nz; 3155 for (i=0; i<n; i++) { 3156 aij->ilen[i] = aij->imax[i]; 3157 } 3158 3159 PetscFunctionReturn(0); 3160 } 3161 EXTERN_C_END 3162 3163 #undef __FUNCT__ 3164 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3165 /*@ 3166 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3167 in the matrix. 3168 3169 Input Parameters: 3170 + mat - the SeqAIJ matrix 3171 - indices - the column indices 3172 3173 Level: advanced 3174 3175 Notes: 3176 This can be called if you have precomputed the nonzero structure of the 3177 matrix and want to provide it to the matrix object to improve the performance 3178 of the MatSetValues() operation. 3179 3180 You MUST have set the correct numbers of nonzeros per row in the call to 3181 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3182 3183 MUST be called before any calls to MatSetValues(); 3184 3185 The indices should start with zero, not one. 3186 3187 @*/ 3188 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3189 { 3190 PetscErrorCode ierr; 3191 3192 PetscFunctionBegin; 3193 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3194 PetscValidPointer(indices,2); 3195 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 3196 PetscFunctionReturn(0); 3197 } 3198 3199 /* ----------------------------------------------------------------------------------------*/ 3200 3201 EXTERN_C_BEGIN 3202 #undef __FUNCT__ 3203 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3204 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3205 { 3206 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3207 PetscErrorCode ierr; 3208 size_t nz = aij->i[mat->rmap->n]; 3209 3210 PetscFunctionBegin; 3211 if (aij->nonew != 1) { 3212 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3213 } 3214 3215 /* allocate space for values if not already there */ 3216 if (!aij->saved_values) { 3217 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3218 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3219 } 3220 3221 /* copy values over */ 3222 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3223 PetscFunctionReturn(0); 3224 } 3225 EXTERN_C_END 3226 3227 #undef __FUNCT__ 3228 #define __FUNCT__ "MatStoreValues" 3229 /*@ 3230 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3231 example, reuse of the linear part of a Jacobian, while recomputing the 3232 nonlinear portion. 3233 3234 Collect on Mat 3235 3236 Input Parameters: 3237 . mat - the matrix (currently only AIJ matrices support this option) 3238 3239 Level: advanced 3240 3241 Common Usage, with SNESSolve(): 3242 $ Create Jacobian matrix 3243 $ Set linear terms into matrix 3244 $ Apply boundary conditions to matrix, at this time matrix must have 3245 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3246 $ boundary conditions again will not change the nonzero structure 3247 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3248 $ ierr = MatStoreValues(mat); 3249 $ Call SNESSetJacobian() with matrix 3250 $ In your Jacobian routine 3251 $ ierr = MatRetrieveValues(mat); 3252 $ Set nonlinear terms in matrix 3253 3254 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3255 $ // build linear portion of Jacobian 3256 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3257 $ ierr = MatStoreValues(mat); 3258 $ loop over nonlinear iterations 3259 $ ierr = MatRetrieveValues(mat); 3260 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3261 $ // call MatAssemblyBegin/End() on matrix 3262 $ Solve linear system with Jacobian 3263 $ endloop 3264 3265 Notes: 3266 Matrix must already be assemblied before calling this routine 3267 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3268 calling this routine. 3269 3270 When this is called multiple times it overwrites the previous set of stored values 3271 and does not allocated additional space. 3272 3273 .seealso: MatRetrieveValues() 3274 3275 @*/ 3276 PetscErrorCode MatStoreValues(Mat mat) 3277 { 3278 PetscErrorCode ierr; 3279 3280 PetscFunctionBegin; 3281 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3282 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3283 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3284 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3285 PetscFunctionReturn(0); 3286 } 3287 3288 EXTERN_C_BEGIN 3289 #undef __FUNCT__ 3290 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3291 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3292 { 3293 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3294 PetscErrorCode ierr; 3295 PetscInt nz = aij->i[mat->rmap->n]; 3296 3297 PetscFunctionBegin; 3298 if (aij->nonew != 1) { 3299 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3300 } 3301 if (!aij->saved_values) { 3302 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3303 } 3304 /* copy values over */ 3305 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3306 PetscFunctionReturn(0); 3307 } 3308 EXTERN_C_END 3309 3310 #undef __FUNCT__ 3311 #define __FUNCT__ "MatRetrieveValues" 3312 /*@ 3313 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3314 example, reuse of the linear part of a Jacobian, while recomputing the 3315 nonlinear portion. 3316 3317 Collect on Mat 3318 3319 Input Parameters: 3320 . mat - the matrix (currently on AIJ matrices support this option) 3321 3322 Level: advanced 3323 3324 .seealso: MatStoreValues() 3325 3326 @*/ 3327 PetscErrorCode MatRetrieveValues(Mat mat) 3328 { 3329 PetscErrorCode ierr; 3330 3331 PetscFunctionBegin; 3332 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3333 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3334 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3335 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3336 PetscFunctionReturn(0); 3337 } 3338 3339 3340 /* --------------------------------------------------------------------------------*/ 3341 #undef __FUNCT__ 3342 #define __FUNCT__ "MatCreateSeqAIJ" 3343 /*@C 3344 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3345 (the default parallel PETSc format). For good matrix assembly performance 3346 the user should preallocate the matrix storage by setting the parameter nz 3347 (or the array nnz). By setting these parameters accurately, performance 3348 during matrix assembly can be increased by more than a factor of 50. 3349 3350 Collective on MPI_Comm 3351 3352 Input Parameters: 3353 + comm - MPI communicator, set to PETSC_COMM_SELF 3354 . m - number of rows 3355 . n - number of columns 3356 . nz - number of nonzeros per row (same for all rows) 3357 - nnz - array containing the number of nonzeros in the various rows 3358 (possibly different for each row) or PETSC_NULL 3359 3360 Output Parameter: 3361 . A - the matrix 3362 3363 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3364 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3365 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3366 3367 Notes: 3368 If nnz is given then nz is ignored 3369 3370 The AIJ format (also called the Yale sparse matrix format or 3371 compressed row storage), is fully compatible with standard Fortran 77 3372 storage. That is, the stored row and column indices can begin at 3373 either one (as in Fortran) or zero. See the users' manual for details. 3374 3375 Specify the preallocated storage with either nz or nnz (not both). 3376 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3377 allocation. For large problems you MUST preallocate memory or you 3378 will get TERRIBLE performance, see the users' manual chapter on matrices. 3379 3380 By default, this format uses inodes (identical nodes) when possible, to 3381 improve numerical efficiency of matrix-vector products and solves. We 3382 search for consecutive rows with the same nonzero structure, thereby 3383 reusing matrix information to achieve increased efficiency. 3384 3385 Options Database Keys: 3386 + -mat_no_inode - Do not use inodes 3387 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3388 3389 Level: intermediate 3390 3391 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3392 3393 @*/ 3394 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3395 { 3396 PetscErrorCode ierr; 3397 3398 PetscFunctionBegin; 3399 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3400 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3401 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3402 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3403 PetscFunctionReturn(0); 3404 } 3405 3406 #undef __FUNCT__ 3407 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3408 /*@C 3409 MatSeqAIJSetPreallocation - For good matrix assembly performance 3410 the user should preallocate the matrix storage by setting the parameter nz 3411 (or the array nnz). By setting these parameters accurately, performance 3412 during matrix assembly can be increased by more than a factor of 50. 3413 3414 Collective on MPI_Comm 3415 3416 Input Parameters: 3417 + B - The matrix-free 3418 . nz - number of nonzeros per row (same for all rows) 3419 - nnz - array containing the number of nonzeros in the various rows 3420 (possibly different for each row) or PETSC_NULL 3421 3422 Notes: 3423 If nnz is given then nz is ignored 3424 3425 The AIJ format (also called the Yale sparse matrix format or 3426 compressed row storage), is fully compatible with standard Fortran 77 3427 storage. That is, the stored row and column indices can begin at 3428 either one (as in Fortran) or zero. See the users' manual for details. 3429 3430 Specify the preallocated storage with either nz or nnz (not both). 3431 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3432 allocation. For large problems you MUST preallocate memory or you 3433 will get TERRIBLE performance, see the users' manual chapter on matrices. 3434 3435 You can call MatGetInfo() to get information on how effective the preallocation was; 3436 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3437 You can also run with the option -info and look for messages with the string 3438 malloc in them to see if additional memory allocation was needed. 3439 3440 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3441 entries or columns indices 3442 3443 By default, this format uses inodes (identical nodes) when possible, to 3444 improve numerical efficiency of matrix-vector products and solves. We 3445 search for consecutive rows with the same nonzero structure, thereby 3446 reusing matrix information to achieve increased efficiency. 3447 3448 Options Database Keys: 3449 + -mat_no_inode - Do not use inodes 3450 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3451 - -mat_aij_oneindex - Internally use indexing starting at 1 3452 rather than 0. Note that when calling MatSetValues(), 3453 the user still MUST index entries starting at 0! 3454 3455 Level: intermediate 3456 3457 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3458 3459 @*/ 3460 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3461 { 3462 PetscErrorCode ierr; 3463 3464 PetscFunctionBegin; 3465 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3466 PetscFunctionReturn(0); 3467 } 3468 3469 EXTERN_C_BEGIN 3470 #undef __FUNCT__ 3471 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3472 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3473 { 3474 Mat_SeqAIJ *b; 3475 PetscBool skipallocation = PETSC_FALSE; 3476 PetscErrorCode ierr; 3477 PetscInt i; 3478 3479 PetscFunctionBegin; 3480 3481 if (nz == MAT_SKIP_ALLOCATION) { 3482 skipallocation = PETSC_TRUE; 3483 nz = 0; 3484 } 3485 3486 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3487 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3488 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3489 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3490 3491 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3492 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3493 if (nnz) { 3494 for (i=0; i<B->rmap->n; i++) { 3495 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]); 3496 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); 3497 } 3498 } 3499 3500 B->preallocated = PETSC_TRUE; 3501 b = (Mat_SeqAIJ*)B->data; 3502 3503 if (!skipallocation) { 3504 if (!b->imax) { 3505 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3506 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3507 } 3508 if (!nnz) { 3509 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3510 else if (nz < 0) nz = 1; 3511 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3512 nz = nz*B->rmap->n; 3513 } else { 3514 nz = 0; 3515 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3516 } 3517 /* b->ilen will count nonzeros in each row so far. */ 3518 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3519 3520 /* allocate the matrix space */ 3521 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3522 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3523 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3524 b->i[0] = 0; 3525 for (i=1; i<B->rmap->n+1; i++) { 3526 b->i[i] = b->i[i-1] + b->imax[i-1]; 3527 } 3528 b->singlemalloc = PETSC_TRUE; 3529 b->free_a = PETSC_TRUE; 3530 b->free_ij = PETSC_TRUE; 3531 } else { 3532 b->free_a = PETSC_FALSE; 3533 b->free_ij = PETSC_FALSE; 3534 } 3535 3536 b->nz = 0; 3537 b->maxnz = nz; 3538 B->info.nz_unneeded = (double)b->maxnz; 3539 PetscFunctionReturn(0); 3540 } 3541 EXTERN_C_END 3542 3543 #undef __FUNCT__ 3544 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3545 /*@ 3546 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3547 3548 Input Parameters: 3549 + B - the matrix 3550 . i - the indices into j for the start of each row (starts with zero) 3551 . j - the column indices for each row (starts with zero) these must be sorted for each row 3552 - v - optional values in the matrix 3553 3554 Level: developer 3555 3556 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3557 3558 .keywords: matrix, aij, compressed row, sparse, sequential 3559 3560 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3561 @*/ 3562 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3563 { 3564 PetscErrorCode ierr; 3565 3566 PetscFunctionBegin; 3567 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3568 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3569 PetscFunctionReturn(0); 3570 } 3571 3572 EXTERN_C_BEGIN 3573 #undef __FUNCT__ 3574 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3575 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3576 { 3577 PetscInt i; 3578 PetscInt m,n; 3579 PetscInt nz; 3580 PetscInt *nnz, nz_max = 0; 3581 PetscScalar *values; 3582 PetscErrorCode ierr; 3583 3584 PetscFunctionBegin; 3585 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3586 3587 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3588 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3589 for(i = 0; i < m; i++) { 3590 nz = Ii[i+1]- Ii[i]; 3591 nz_max = PetscMax(nz_max, nz); 3592 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3593 nnz[i] = nz; 3594 } 3595 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3596 ierr = PetscFree(nnz);CHKERRQ(ierr); 3597 3598 if (v) { 3599 values = (PetscScalar*) v; 3600 } else { 3601 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3602 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3603 } 3604 3605 for(i = 0; i < m; i++) { 3606 nz = Ii[i+1] - Ii[i]; 3607 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3608 } 3609 3610 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3611 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3612 3613 if (!v) { 3614 ierr = PetscFree(values);CHKERRQ(ierr); 3615 } 3616 PetscFunctionReturn(0); 3617 } 3618 EXTERN_C_END 3619 3620 #include <../src/mat/impls/dense/seq/dense.h> 3621 #include <private/petscaxpy.h> 3622 3623 #undef __FUNCT__ 3624 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3625 /* 3626 Computes (B'*A')' since computing B*A directly is untenable 3627 3628 n p p 3629 ( ) ( ) ( ) 3630 m ( A ) * n ( B ) = m ( C ) 3631 ( ) ( ) ( ) 3632 3633 */ 3634 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3635 { 3636 PetscErrorCode ierr; 3637 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3638 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3639 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3640 PetscInt i,n,m,q,p; 3641 const PetscInt *ii,*idx; 3642 const PetscScalar *b,*a,*a_q; 3643 PetscScalar *c,*c_q; 3644 3645 PetscFunctionBegin; 3646 m = A->rmap->n; 3647 n = A->cmap->n; 3648 p = B->cmap->n; 3649 a = sub_a->v; 3650 b = sub_b->a; 3651 c = sub_c->v; 3652 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3653 3654 ii = sub_b->i; 3655 idx = sub_b->j; 3656 for (i=0; i<n; i++) { 3657 q = ii[i+1] - ii[i]; 3658 while (q-->0) { 3659 c_q = c + m*(*idx); 3660 a_q = a + m*i; 3661 PetscAXPY(c_q,*b,a_q,m); 3662 idx++; 3663 b++; 3664 } 3665 } 3666 PetscFunctionReturn(0); 3667 } 3668 3669 #undef __FUNCT__ 3670 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3671 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3672 { 3673 PetscErrorCode ierr; 3674 PetscInt m=A->rmap->n,n=B->cmap->n; 3675 Mat Cmat; 3676 3677 PetscFunctionBegin; 3678 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); 3679 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3680 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3681 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3682 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3683 Cmat->assembled = PETSC_TRUE; 3684 *C = Cmat; 3685 PetscFunctionReturn(0); 3686 } 3687 3688 /* ----------------------------------------------------------------*/ 3689 #undef __FUNCT__ 3690 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3691 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3692 { 3693 PetscErrorCode ierr; 3694 3695 PetscFunctionBegin; 3696 if (scall == MAT_INITIAL_MATRIX){ 3697 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3698 } 3699 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3700 PetscFunctionReturn(0); 3701 } 3702 3703 3704 /*MC 3705 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3706 based on compressed sparse row format. 3707 3708 Options Database Keys: 3709 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3710 3711 Level: beginner 3712 3713 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3714 M*/ 3715 3716 EXTERN_C_BEGIN 3717 #if defined(PETSC_HAVE_PASTIX) 3718 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3719 #endif 3720 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3721 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3722 #endif 3723 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3724 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3725 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3726 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3727 #if defined(PETSC_HAVE_MUMPS) 3728 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3729 #endif 3730 #if defined(PETSC_HAVE_SUPERLU) 3731 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3732 #endif 3733 #if defined(PETSC_HAVE_SUPERLU_DIST) 3734 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3735 #endif 3736 #if defined(PETSC_HAVE_SPOOLES) 3737 extern PetscErrorCode MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3738 #endif 3739 #if defined(PETSC_HAVE_UMFPACK) 3740 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3741 #endif 3742 #if defined(PETSC_HAVE_CHOLMOD) 3743 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3744 #endif 3745 #if defined(PETSC_HAVE_LUSOL) 3746 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3747 #endif 3748 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3749 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3750 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3751 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3752 #endif 3753 EXTERN_C_END 3754 3755 EXTERN_C_BEGIN 3756 #undef __FUNCT__ 3757 #define __FUNCT__ "MatCreate_SeqAIJ" 3758 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3759 { 3760 Mat_SeqAIJ *b; 3761 PetscErrorCode ierr; 3762 PetscMPIInt size; 3763 3764 PetscFunctionBegin; 3765 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3766 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3767 3768 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3769 B->data = (void*)b; 3770 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3771 b->row = 0; 3772 b->col = 0; 3773 b->icol = 0; 3774 b->reallocs = 0; 3775 b->ignorezeroentries = PETSC_FALSE; 3776 b->roworiented = PETSC_TRUE; 3777 b->nonew = 0; 3778 b->diag = 0; 3779 b->solve_work = 0; 3780 B->spptr = 0; 3781 b->saved_values = 0; 3782 b->idiag = 0; 3783 b->mdiag = 0; 3784 b->ssor_work = 0; 3785 b->omega = 1.0; 3786 b->fshift = 0.0; 3787 b->idiagvalid = PETSC_FALSE; 3788 b->ibdiagvalid = PETSC_FALSE; 3789 b->keepnonzeropattern = PETSC_FALSE; 3790 b->xtoy = 0; 3791 b->XtoY = 0; 3792 B->same_nonzero = PETSC_FALSE; 3793 3794 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3795 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3796 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3797 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3798 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3799 #endif 3800 #if defined(PETSC_HAVE_PASTIX) 3801 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3802 #endif 3803 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3804 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3805 #endif 3806 #if defined(PETSC_HAVE_SUPERLU) 3807 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3808 #endif 3809 #if defined(PETSC_HAVE_SUPERLU_DIST) 3810 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3811 #endif 3812 #if defined(PETSC_HAVE_SPOOLES) 3813 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C","MatGetFactor_seqaij_spooles",MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3814 #endif 3815 #if defined(PETSC_HAVE_MUMPS) 3816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3817 #endif 3818 #if defined(PETSC_HAVE_UMFPACK) 3819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3820 #endif 3821 #if defined(PETSC_HAVE_CHOLMOD) 3822 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3823 #endif 3824 #if defined(PETSC_HAVE_LUSOL) 3825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_seqaij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3826 #endif 3827 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3828 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3829 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3830 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3831 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3832 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3833 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3834 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3835 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3836 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3837 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3838 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3839 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3840 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3841 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3842 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3843 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3844 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3845 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3846 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3847 PetscFunctionReturn(0); 3848 } 3849 EXTERN_C_END 3850 3851 #if defined(PETSC_HAVE_PTHREADCLASSES) 3852 EXTERN_C_BEGIN 3853 #undef __FUNCT__ 3854 #define __FUNCT__ "MatCreate_SeqPThreadAIJ" 3855 PetscErrorCode MatCreate_SeqPThreadAIJ(Mat B) 3856 { 3857 PetscErrorCode ierr; 3858 3859 PetscFunctionBegin; 3860 ierr = MatCreate_SeqAIJ(B); 3861 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3862 B->ops->mult = MatMult_SeqPThreadAIJ; 3863 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQPTHREADAIJ);CHKERRQ(ierr); 3864 PetscFunctionReturn(0); 3865 } 3866 EXTERN_C_END 3867 #endif 3868 3869 #undef __FUNCT__ 3870 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3871 /* 3872 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3873 */ 3874 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3875 { 3876 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3877 PetscErrorCode ierr; 3878 PetscInt i,m = A->rmap->n; 3879 3880 PetscFunctionBegin; 3881 c = (Mat_SeqAIJ*)C->data; 3882 3883 C->factortype = A->factortype; 3884 c->row = 0; 3885 c->col = 0; 3886 c->icol = 0; 3887 c->reallocs = 0; 3888 3889 C->assembled = PETSC_TRUE; 3890 3891 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3892 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3893 3894 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3895 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3896 for (i=0; i<m; i++) { 3897 c->imax[i] = a->imax[i]; 3898 c->ilen[i] = a->ilen[i]; 3899 } 3900 3901 /* allocate the matrix space */ 3902 if (mallocmatspace){ 3903 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3904 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3905 c->singlemalloc = PETSC_TRUE; 3906 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3907 if (m > 0) { 3908 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3909 if (cpvalues == MAT_COPY_VALUES) { 3910 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3911 } else { 3912 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3913 } 3914 } 3915 } 3916 3917 c->ignorezeroentries = a->ignorezeroentries; 3918 c->roworiented = a->roworiented; 3919 c->nonew = a->nonew; 3920 if (a->diag) { 3921 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3922 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3923 for (i=0; i<m; i++) { 3924 c->diag[i] = a->diag[i]; 3925 } 3926 } else c->diag = 0; 3927 c->solve_work = 0; 3928 c->saved_values = 0; 3929 c->idiag = 0; 3930 c->ssor_work = 0; 3931 c->keepnonzeropattern = a->keepnonzeropattern; 3932 c->free_a = PETSC_TRUE; 3933 c->free_ij = PETSC_TRUE; 3934 c->xtoy = 0; 3935 c->XtoY = 0; 3936 3937 c->nz = a->nz; 3938 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3939 C->preallocated = PETSC_TRUE; 3940 3941 c->compressedrow.use = a->compressedrow.use; 3942 c->compressedrow.nrows = a->compressedrow.nrows; 3943 c->compressedrow.check = a->compressedrow.check; 3944 if (a->compressedrow.use){ 3945 i = a->compressedrow.nrows; 3946 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3947 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3948 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3949 } else { 3950 c->compressedrow.use = PETSC_FALSE; 3951 c->compressedrow.i = PETSC_NULL; 3952 c->compressedrow.rindex = PETSC_NULL; 3953 } 3954 C->same_nonzero = A->same_nonzero; 3955 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3956 3957 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3958 PetscFunctionReturn(0); 3959 } 3960 3961 #undef __FUNCT__ 3962 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3963 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3964 { 3965 PetscErrorCode ierr; 3966 3967 PetscFunctionBegin; 3968 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3969 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3970 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3971 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3972 PetscFunctionReturn(0); 3973 } 3974 3975 #undef __FUNCT__ 3976 #define __FUNCT__ "MatLoad_SeqAIJ" 3977 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3978 { 3979 Mat_SeqAIJ *a; 3980 PetscErrorCode ierr; 3981 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3982 int fd; 3983 PetscMPIInt size; 3984 MPI_Comm comm; 3985 PetscInt bs = 1; 3986 3987 PetscFunctionBegin; 3988 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3989 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3990 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3991 3992 ierr = PetscOptionsBegin(comm,PETSC_NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 3993 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 3994 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3995 3996 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3997 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3998 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3999 M = header[1]; N = header[2]; nz = header[3]; 4000 4001 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4002 4003 /* read in row lengths */ 4004 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4005 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4006 4007 /* check if sum of rowlengths is same as nz */ 4008 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4009 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); 4010 4011 /* set global size if not set already*/ 4012 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4013 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4014 } else { 4015 /* if sizes and type are already set, check if the vector global sizes are correct */ 4016 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4017 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); 4018 } 4019 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4020 a = (Mat_SeqAIJ*)newMat->data; 4021 4022 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4023 4024 /* read in nonzero values */ 4025 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4026 4027 /* set matrix "i" values */ 4028 a->i[0] = 0; 4029 for (i=1; i<= M; i++) { 4030 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4031 a->ilen[i-1] = rowlengths[i-1]; 4032 } 4033 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4034 4035 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4036 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4037 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4038 PetscFunctionReturn(0); 4039 } 4040 4041 #undef __FUNCT__ 4042 #define __FUNCT__ "MatEqual_SeqAIJ" 4043 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4044 { 4045 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 4046 PetscErrorCode ierr; 4047 #if defined(PETSC_USE_COMPLEX) 4048 PetscInt k; 4049 #endif 4050 4051 PetscFunctionBegin; 4052 /* If the matrix dimensions are not equal,or no of nonzeros */ 4053 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4054 *flg = PETSC_FALSE; 4055 PetscFunctionReturn(0); 4056 } 4057 4058 /* if the a->i are the same */ 4059 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4060 if (!*flg) PetscFunctionReturn(0); 4061 4062 /* if a->j are the same */ 4063 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4064 if (!*flg) PetscFunctionReturn(0); 4065 4066 /* if a->a are the same */ 4067 #if defined(PETSC_USE_COMPLEX) 4068 for (k=0; k<a->nz; k++){ 4069 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 4070 *flg = PETSC_FALSE; 4071 PetscFunctionReturn(0); 4072 } 4073 } 4074 #else 4075 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4076 #endif 4077 PetscFunctionReturn(0); 4078 } 4079 4080 #undef __FUNCT__ 4081 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4082 /*@ 4083 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4084 provided by the user. 4085 4086 Collective on MPI_Comm 4087 4088 Input Parameters: 4089 + comm - must be an MPI communicator of size 1 4090 . m - number of rows 4091 . n - number of columns 4092 . i - row indices 4093 . j - column indices 4094 - a - matrix values 4095 4096 Output Parameter: 4097 . mat - the matrix 4098 4099 Level: intermediate 4100 4101 Notes: 4102 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4103 once the matrix is destroyed and not before 4104 4105 You cannot set new nonzero locations into this matrix, that will generate an error. 4106 4107 The i and j indices are 0 based 4108 4109 The format which is used for the sparse matrix input, is equivalent to a 4110 row-major ordering.. i.e for the following matrix, the input data expected is 4111 as shown: 4112 4113 1 0 0 4114 2 0 3 4115 4 5 6 4116 4117 i = {0,1,3,6} [size = nrow+1 = 3+1] 4118 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4119 v = {1,2,3,4,5,6} [size = nz = 6] 4120 4121 4122 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4123 4124 @*/ 4125 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 4126 { 4127 PetscErrorCode ierr; 4128 PetscInt ii; 4129 Mat_SeqAIJ *aij; 4130 #if defined(PETSC_USE_DEBUG) 4131 PetscInt jj; 4132 #endif 4133 4134 PetscFunctionBegin; 4135 if (i[0]) { 4136 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4137 } 4138 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4139 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4140 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4141 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4142 aij = (Mat_SeqAIJ*)(*mat)->data; 4143 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4144 4145 aij->i = i; 4146 aij->j = j; 4147 aij->a = a; 4148 aij->singlemalloc = PETSC_FALSE; 4149 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4150 aij->free_a = PETSC_FALSE; 4151 aij->free_ij = PETSC_FALSE; 4152 4153 for (ii=0; ii<m; ii++) { 4154 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4155 #if defined(PETSC_USE_DEBUG) 4156 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]); 4157 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4158 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); 4159 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 4160 } 4161 #endif 4162 } 4163 #if defined(PETSC_USE_DEBUG) 4164 for (ii=0; ii<aij->i[m]; ii++) { 4165 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4166 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]); 4167 } 4168 #endif 4169 4170 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4171 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4172 PetscFunctionReturn(0); 4173 } 4174 4175 #undef __FUNCT__ 4176 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4177 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4178 { 4179 PetscErrorCode ierr; 4180 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4181 4182 PetscFunctionBegin; 4183 if (coloring->ctype == IS_COLORING_GLOBAL) { 4184 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4185 a->coloring = coloring; 4186 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4187 PetscInt i,*larray; 4188 ISColoring ocoloring; 4189 ISColoringValue *colors; 4190 4191 /* set coloring for diagonal portion */ 4192 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4193 for (i=0; i<A->cmap->n; i++) { 4194 larray[i] = i; 4195 } 4196 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 4197 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4198 for (i=0; i<A->cmap->n; i++) { 4199 colors[i] = coloring->colors[larray[i]]; 4200 } 4201 ierr = PetscFree(larray);CHKERRQ(ierr); 4202 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4203 a->coloring = ocoloring; 4204 } 4205 PetscFunctionReturn(0); 4206 } 4207 4208 #if defined(PETSC_HAVE_ADIC) 4209 EXTERN_C_BEGIN 4210 #include <adic/ad_utils.h> 4211 EXTERN_C_END 4212 4213 #undef __FUNCT__ 4214 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 4215 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 4216 { 4217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4218 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 4219 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 4220 ISColoringValue *color; 4221 4222 PetscFunctionBegin; 4223 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4224 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 4225 color = a->coloring->colors; 4226 /* loop over rows */ 4227 for (i=0; i<m; i++) { 4228 nz = ii[i+1] - ii[i]; 4229 /* loop over columns putting computed value into matrix */ 4230 for (j=0; j<nz; j++) { 4231 *v++ = values[color[*jj++]]; 4232 } 4233 values += nlen; /* jump to next row of derivatives */ 4234 } 4235 PetscFunctionReturn(0); 4236 } 4237 #endif 4238 4239 #undef __FUNCT__ 4240 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4241 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4242 { 4243 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4244 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4245 MatScalar *v = a->a; 4246 PetscScalar *values = (PetscScalar *)advalues; 4247 ISColoringValue *color; 4248 4249 PetscFunctionBegin; 4250 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4251 color = a->coloring->colors; 4252 /* loop over rows */ 4253 for (i=0; i<m; i++) { 4254 nz = ii[i+1] - ii[i]; 4255 /* loop over columns putting computed value into matrix */ 4256 for (j=0; j<nz; j++) { 4257 *v++ = values[color[*jj++]]; 4258 } 4259 values += nl; /* jump to next row of derivatives */ 4260 } 4261 PetscFunctionReturn(0); 4262 } 4263 4264 /* 4265 Special version for direct calls from Fortran 4266 */ 4267 #include <private/fortranimpl.h> 4268 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4269 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4270 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4271 #define matsetvaluesseqaij_ matsetvaluesseqaij 4272 #endif 4273 4274 /* Change these macros so can be used in void function */ 4275 #undef CHKERRQ 4276 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 4277 #undef SETERRQ2 4278 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4279 4280 EXTERN_C_BEGIN 4281 #undef __FUNCT__ 4282 #define __FUNCT__ "matsetvaluesseqaij_" 4283 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4284 { 4285 Mat A = *AA; 4286 PetscInt m = *mm, n = *nn; 4287 InsertMode is = *isis; 4288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4289 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4290 PetscInt *imax,*ai,*ailen; 4291 PetscErrorCode ierr; 4292 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4293 MatScalar *ap,value,*aa; 4294 PetscBool ignorezeroentries = a->ignorezeroentries; 4295 PetscBool roworiented = a->roworiented; 4296 4297 PetscFunctionBegin; 4298 ierr = MatPreallocated(A);CHKERRQ(ierr); 4299 imax = a->imax; 4300 ai = a->i; 4301 ailen = a->ilen; 4302 aj = a->j; 4303 aa = a->a; 4304 4305 for (k=0; k<m; k++) { /* loop over added rows */ 4306 row = im[k]; 4307 if (row < 0) continue; 4308 #if defined(PETSC_USE_DEBUG) 4309 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4310 #endif 4311 rp = aj + ai[row]; ap = aa + ai[row]; 4312 rmax = imax[row]; nrow = ailen[row]; 4313 low = 0; 4314 high = nrow; 4315 for (l=0; l<n; l++) { /* loop over added columns */ 4316 if (in[l] < 0) continue; 4317 #if defined(PETSC_USE_DEBUG) 4318 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4319 #endif 4320 col = in[l]; 4321 if (roworiented) { 4322 value = v[l + k*n]; 4323 } else { 4324 value = v[k + l*m]; 4325 } 4326 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4327 4328 if (col <= lastcol) low = 0; else high = nrow; 4329 lastcol = col; 4330 while (high-low > 5) { 4331 t = (low+high)/2; 4332 if (rp[t] > col) high = t; 4333 else low = t; 4334 } 4335 for (i=low; i<high; i++) { 4336 if (rp[i] > col) break; 4337 if (rp[i] == col) { 4338 if (is == ADD_VALUES) ap[i] += value; 4339 else ap[i] = value; 4340 goto noinsert; 4341 } 4342 } 4343 if (value == 0.0 && ignorezeroentries) goto noinsert; 4344 if (nonew == 1) goto noinsert; 4345 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4346 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4347 N = nrow++ - 1; a->nz++; high++; 4348 /* shift up all the later entries in this row */ 4349 for (ii=N; ii>=i; ii--) { 4350 rp[ii+1] = rp[ii]; 4351 ap[ii+1] = ap[ii]; 4352 } 4353 rp[i] = col; 4354 ap[i] = value; 4355 noinsert:; 4356 low = i + 1; 4357 } 4358 ailen[row] = nrow; 4359 } 4360 A->same_nonzero = PETSC_FALSE; 4361 PetscFunctionReturnVoid(); 4362 } 4363 EXTERN_C_END 4364 4365