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