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