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 #include <../src/mat/blocktranspose.h> 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 ierr = Kernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2964 diag += 4; 2965 } 2966 break; 2967 case 3: 2968 for (i=0; i<mbs; i++) { 2969 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2970 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2971 ierr = Kernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 2972 ierr = Kernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2973 diag += 9; 2974 } 2975 break; 2976 case 4: 2977 for (i=0; i<mbs; i++) { 2978 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2979 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2980 ierr = Kernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 2981 ierr = Kernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2982 diag += 16; 2983 } 2984 break; 2985 case 5: 2986 for (i=0; i<mbs; i++) { 2987 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2988 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2989 ierr = Kernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 2990 ierr = Kernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 2991 diag += 25; 2992 } 2993 break; 2994 case 6: 2995 for (i=0; i<mbs; i++) { 2996 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; 2997 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2998 ierr = Kernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 2999 ierr = Kernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3000 diag += 36; 3001 } 3002 break; 3003 case 7: 3004 for (i=0; i<mbs; i++) { 3005 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; 3006 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3007 ierr = Kernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3008 ierr = Kernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3009 diag += 49; 3010 } 3011 break; 3012 default: 3013 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 3014 for (i=0; i<mbs; i++) { 3015 for (j=0; j<bs; j++) { 3016 IJ[j] = bs*i + j; 3017 } 3018 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3019 ierr = Kernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3020 ierr = Kernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3021 diag += bs2; 3022 } 3023 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3024 } 3025 a->ibdiagvalid = PETSC_TRUE; 3026 PetscFunctionReturn(0); 3027 } 3028 3029 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 3030 /* -------------------------------------------------------------------*/ 3031 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3032 MatGetRow_SeqAIJ, 3033 MatRestoreRow_SeqAIJ, 3034 MatMult_SeqAIJ, 3035 /* 4*/ MatMultAdd_SeqAIJ, 3036 MatMultTranspose_SeqAIJ, 3037 MatMultTransposeAdd_SeqAIJ, 3038 0, 3039 0, 3040 0, 3041 /*10*/ 0, 3042 MatLUFactor_SeqAIJ, 3043 0, 3044 MatSOR_SeqAIJ, 3045 MatTranspose_SeqAIJ, 3046 /*15*/ MatGetInfo_SeqAIJ, 3047 MatEqual_SeqAIJ, 3048 MatGetDiagonal_SeqAIJ, 3049 MatDiagonalScale_SeqAIJ, 3050 MatNorm_SeqAIJ, 3051 /*20*/ 0, 3052 MatAssemblyEnd_SeqAIJ, 3053 MatSetOption_SeqAIJ, 3054 MatZeroEntries_SeqAIJ, 3055 /*24*/ MatZeroRows_SeqAIJ, 3056 0, 3057 0, 3058 0, 3059 0, 3060 /*29*/ MatSetUpPreallocation_SeqAIJ, 3061 0, 3062 0, 3063 MatGetArray_SeqAIJ, 3064 MatRestoreArray_SeqAIJ, 3065 /*34*/ MatDuplicate_SeqAIJ, 3066 0, 3067 0, 3068 MatILUFactor_SeqAIJ, 3069 0, 3070 /*39*/ MatAXPY_SeqAIJ, 3071 MatGetSubMatrices_SeqAIJ, 3072 MatIncreaseOverlap_SeqAIJ, 3073 MatGetValues_SeqAIJ, 3074 MatCopy_SeqAIJ, 3075 /*44*/ MatGetRowMax_SeqAIJ, 3076 MatScale_SeqAIJ, 3077 0, 3078 MatDiagonalSet_SeqAIJ, 3079 MatZeroRowsColumns_SeqAIJ, 3080 /*49*/ MatSetBlockSize_SeqAIJ, 3081 MatGetRowIJ_SeqAIJ, 3082 MatRestoreRowIJ_SeqAIJ, 3083 MatGetColumnIJ_SeqAIJ, 3084 MatRestoreColumnIJ_SeqAIJ, 3085 /*54*/ MatFDColoringCreate_SeqAIJ, 3086 0, 3087 0, 3088 MatPermute_SeqAIJ, 3089 0, 3090 /*59*/ 0, 3091 MatDestroy_SeqAIJ, 3092 MatView_SeqAIJ, 3093 0, 3094 0, 3095 /*64*/ 0, 3096 0, 3097 0, 3098 0, 3099 0, 3100 /*69*/ MatGetRowMaxAbs_SeqAIJ, 3101 MatGetRowMinAbs_SeqAIJ, 3102 0, 3103 MatSetColoring_SeqAIJ, 3104 #if defined(PETSC_HAVE_ADIC) 3105 MatSetValuesAdic_SeqAIJ, 3106 #else 3107 0, 3108 #endif 3109 /*74*/ MatSetValuesAdifor_SeqAIJ, 3110 MatFDColoringApply_AIJ, 3111 0, 3112 0, 3113 0, 3114 /*79*/ MatFindZeroDiagonals_SeqAIJ, 3115 0, 3116 0, 3117 0, 3118 MatLoad_SeqAIJ, 3119 /*84*/ MatIsSymmetric_SeqAIJ, 3120 MatIsHermitian_SeqAIJ, 3121 0, 3122 0, 3123 0, 3124 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 3125 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3126 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3127 MatPtAP_Basic, 3128 MatPtAPSymbolic_SeqAIJ, 3129 /*94*/ MatPtAPNumeric_SeqAIJ, 3130 MatMatMultTranspose_SeqAIJ_SeqAIJ, 3131 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 3132 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 3133 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 3134 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3135 0, 3136 0, 3137 MatConjugate_SeqAIJ, 3138 0, 3139 /*104*/MatSetValuesRow_SeqAIJ, 3140 MatRealPart_SeqAIJ, 3141 MatImaginaryPart_SeqAIJ, 3142 0, 3143 0, 3144 /*109*/MatMatSolve_SeqAIJ, 3145 0, 3146 MatGetRowMin_SeqAIJ, 3147 0, 3148 MatMissingDiagonal_SeqAIJ, 3149 /*114*/0, 3150 0, 3151 0, 3152 0, 3153 0, 3154 /*119*/0, 3155 0, 3156 0, 3157 0, 3158 MatGetMultiProcBlock_SeqAIJ, 3159 /*124*/MatFindNonzeroRows_SeqAIJ, 3160 MatGetColumnNorms_SeqAIJ, 3161 MatInvertBlockDiagonal_SeqAIJ, 3162 0, 3163 0, 3164 /*129*/0 3165 }; 3166 3167 EXTERN_C_BEGIN 3168 #undef __FUNCT__ 3169 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3170 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3171 { 3172 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3173 PetscInt i,nz,n; 3174 3175 PetscFunctionBegin; 3176 3177 nz = aij->maxnz; 3178 n = mat->rmap->n; 3179 for (i=0; i<nz; i++) { 3180 aij->j[i] = indices[i]; 3181 } 3182 aij->nz = nz; 3183 for (i=0; i<n; i++) { 3184 aij->ilen[i] = aij->imax[i]; 3185 } 3186 3187 PetscFunctionReturn(0); 3188 } 3189 EXTERN_C_END 3190 3191 #undef __FUNCT__ 3192 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3193 /*@ 3194 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3195 in the matrix. 3196 3197 Input Parameters: 3198 + mat - the SeqAIJ matrix 3199 - indices - the column indices 3200 3201 Level: advanced 3202 3203 Notes: 3204 This can be called if you have precomputed the nonzero structure of the 3205 matrix and want to provide it to the matrix object to improve the performance 3206 of the MatSetValues() operation. 3207 3208 You MUST have set the correct numbers of nonzeros per row in the call to 3209 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3210 3211 MUST be called before any calls to MatSetValues(); 3212 3213 The indices should start with zero, not one. 3214 3215 @*/ 3216 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3217 { 3218 PetscErrorCode ierr; 3219 3220 PetscFunctionBegin; 3221 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3222 PetscValidPointer(indices,2); 3223 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 3224 PetscFunctionReturn(0); 3225 } 3226 3227 /* ----------------------------------------------------------------------------------------*/ 3228 3229 EXTERN_C_BEGIN 3230 #undef __FUNCT__ 3231 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3232 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3233 { 3234 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3235 PetscErrorCode ierr; 3236 size_t nz = aij->i[mat->rmap->n]; 3237 3238 PetscFunctionBegin; 3239 if (aij->nonew != 1) { 3240 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3241 } 3242 3243 /* allocate space for values if not already there */ 3244 if (!aij->saved_values) { 3245 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3246 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3247 } 3248 3249 /* copy values over */ 3250 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3251 PetscFunctionReturn(0); 3252 } 3253 EXTERN_C_END 3254 3255 #undef __FUNCT__ 3256 #define __FUNCT__ "MatStoreValues" 3257 /*@ 3258 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3259 example, reuse of the linear part of a Jacobian, while recomputing the 3260 nonlinear portion. 3261 3262 Collect on Mat 3263 3264 Input Parameters: 3265 . mat - the matrix (currently only AIJ matrices support this option) 3266 3267 Level: advanced 3268 3269 Common Usage, with SNESSolve(): 3270 $ Create Jacobian matrix 3271 $ Set linear terms into matrix 3272 $ Apply boundary conditions to matrix, at this time matrix must have 3273 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3274 $ boundary conditions again will not change the nonzero structure 3275 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3276 $ ierr = MatStoreValues(mat); 3277 $ Call SNESSetJacobian() with matrix 3278 $ In your Jacobian routine 3279 $ ierr = MatRetrieveValues(mat); 3280 $ Set nonlinear terms in matrix 3281 3282 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3283 $ // build linear portion of Jacobian 3284 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3285 $ ierr = MatStoreValues(mat); 3286 $ loop over nonlinear iterations 3287 $ ierr = MatRetrieveValues(mat); 3288 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3289 $ // call MatAssemblyBegin/End() on matrix 3290 $ Solve linear system with Jacobian 3291 $ endloop 3292 3293 Notes: 3294 Matrix must already be assemblied before calling this routine 3295 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3296 calling this routine. 3297 3298 When this is called multiple times it overwrites the previous set of stored values 3299 and does not allocated additional space. 3300 3301 .seealso: MatRetrieveValues() 3302 3303 @*/ 3304 PetscErrorCode MatStoreValues(Mat mat) 3305 { 3306 PetscErrorCode ierr; 3307 3308 PetscFunctionBegin; 3309 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3310 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3311 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3312 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3313 PetscFunctionReturn(0); 3314 } 3315 3316 EXTERN_C_BEGIN 3317 #undef __FUNCT__ 3318 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3319 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3320 { 3321 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3322 PetscErrorCode ierr; 3323 PetscInt nz = aij->i[mat->rmap->n]; 3324 3325 PetscFunctionBegin; 3326 if (aij->nonew != 1) { 3327 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3328 } 3329 if (!aij->saved_values) { 3330 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3331 } 3332 /* copy values over */ 3333 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3334 PetscFunctionReturn(0); 3335 } 3336 EXTERN_C_END 3337 3338 #undef __FUNCT__ 3339 #define __FUNCT__ "MatRetrieveValues" 3340 /*@ 3341 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3342 example, reuse of the linear part of a Jacobian, while recomputing the 3343 nonlinear portion. 3344 3345 Collect on Mat 3346 3347 Input Parameters: 3348 . mat - the matrix (currently on AIJ matrices support this option) 3349 3350 Level: advanced 3351 3352 .seealso: MatStoreValues() 3353 3354 @*/ 3355 PetscErrorCode MatRetrieveValues(Mat mat) 3356 { 3357 PetscErrorCode ierr; 3358 3359 PetscFunctionBegin; 3360 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3361 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3362 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3363 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3364 PetscFunctionReturn(0); 3365 } 3366 3367 3368 /* --------------------------------------------------------------------------------*/ 3369 #undef __FUNCT__ 3370 #define __FUNCT__ "MatCreateSeqAIJ" 3371 /*@C 3372 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3373 (the default parallel PETSc format). For good matrix assembly performance 3374 the user should preallocate the matrix storage by setting the parameter nz 3375 (or the array nnz). By setting these parameters accurately, performance 3376 during matrix assembly can be increased by more than a factor of 50. 3377 3378 Collective on MPI_Comm 3379 3380 Input Parameters: 3381 + comm - MPI communicator, set to PETSC_COMM_SELF 3382 . m - number of rows 3383 . n - number of columns 3384 . nz - number of nonzeros per row (same for all rows) 3385 - nnz - array containing the number of nonzeros in the various rows 3386 (possibly different for each row) or PETSC_NULL 3387 3388 Output Parameter: 3389 . A - the matrix 3390 3391 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3392 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3393 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3394 3395 Notes: 3396 If nnz is given then nz is ignored 3397 3398 The AIJ format (also called the Yale sparse matrix format or 3399 compressed row storage), is fully compatible with standard Fortran 77 3400 storage. That is, the stored row and column indices can begin at 3401 either one (as in Fortran) or zero. See the users' manual for details. 3402 3403 Specify the preallocated storage with either nz or nnz (not both). 3404 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3405 allocation. For large problems you MUST preallocate memory or you 3406 will get TERRIBLE performance, see the users' manual chapter on matrices. 3407 3408 By default, this format uses inodes (identical nodes) when possible, to 3409 improve numerical efficiency of matrix-vector products and solves. We 3410 search for consecutive rows with the same nonzero structure, thereby 3411 reusing matrix information to achieve increased efficiency. 3412 3413 Options Database Keys: 3414 + -mat_no_inode - Do not use inodes 3415 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3416 3417 Level: intermediate 3418 3419 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3420 3421 @*/ 3422 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3423 { 3424 PetscErrorCode ierr; 3425 3426 PetscFunctionBegin; 3427 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3428 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3429 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3430 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3431 PetscFunctionReturn(0); 3432 } 3433 3434 #undef __FUNCT__ 3435 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3436 /*@C 3437 MatSeqAIJSetPreallocation - For good matrix assembly performance 3438 the user should preallocate the matrix storage by setting the parameter nz 3439 (or the array nnz). By setting these parameters accurately, performance 3440 during matrix assembly can be increased by more than a factor of 50. 3441 3442 Collective on MPI_Comm 3443 3444 Input Parameters: 3445 + B - The matrix-free 3446 . nz - number of nonzeros per row (same for all rows) 3447 - nnz - array containing the number of nonzeros in the various rows 3448 (possibly different for each row) or PETSC_NULL 3449 3450 Notes: 3451 If nnz is given then nz is ignored 3452 3453 The AIJ format (also called the Yale sparse matrix format or 3454 compressed row storage), is fully compatible with standard Fortran 77 3455 storage. That is, the stored row and column indices can begin at 3456 either one (as in Fortran) or zero. See the users' manual for details. 3457 3458 Specify the preallocated storage with either nz or nnz (not both). 3459 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3460 allocation. For large problems you MUST preallocate memory or you 3461 will get TERRIBLE performance, see the users' manual chapter on matrices. 3462 3463 You can call MatGetInfo() to get information on how effective the preallocation was; 3464 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3465 You can also run with the option -info and look for messages with the string 3466 malloc in them to see if additional memory allocation was needed. 3467 3468 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3469 entries or columns indices 3470 3471 By default, this format uses inodes (identical nodes) when possible, to 3472 improve numerical efficiency of matrix-vector products and solves. We 3473 search for consecutive rows with the same nonzero structure, thereby 3474 reusing matrix information to achieve increased efficiency. 3475 3476 Options Database Keys: 3477 + -mat_no_inode - Do not use inodes 3478 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3479 - -mat_aij_oneindex - Internally use indexing starting at 1 3480 rather than 0. Note that when calling MatSetValues(), 3481 the user still MUST index entries starting at 0! 3482 3483 Level: intermediate 3484 3485 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3486 3487 @*/ 3488 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3489 { 3490 PetscErrorCode ierr; 3491 3492 PetscFunctionBegin; 3493 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3494 PetscFunctionReturn(0); 3495 } 3496 3497 EXTERN_C_BEGIN 3498 #undef __FUNCT__ 3499 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3500 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3501 { 3502 Mat_SeqAIJ *b; 3503 PetscBool skipallocation = PETSC_FALSE; 3504 PetscErrorCode ierr; 3505 PetscInt i; 3506 3507 PetscFunctionBegin; 3508 3509 if (nz == MAT_SKIP_ALLOCATION) { 3510 skipallocation = PETSC_TRUE; 3511 nz = 0; 3512 } 3513 3514 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3515 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3516 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3517 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3518 3519 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3520 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3521 if (nnz) { 3522 for (i=0; i<B->rmap->n; i++) { 3523 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]); 3524 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); 3525 } 3526 } 3527 3528 B->preallocated = PETSC_TRUE; 3529 b = (Mat_SeqAIJ*)B->data; 3530 3531 if (!skipallocation) { 3532 if (!b->imax) { 3533 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3534 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3535 } 3536 if (!nnz) { 3537 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3538 else if (nz < 0) nz = 1; 3539 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3540 nz = nz*B->rmap->n; 3541 } else { 3542 nz = 0; 3543 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3544 } 3545 /* b->ilen will count nonzeros in each row so far. */ 3546 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3547 3548 /* allocate the matrix space */ 3549 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3550 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3551 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3552 b->i[0] = 0; 3553 for (i=1; i<B->rmap->n+1; i++) { 3554 b->i[i] = b->i[i-1] + b->imax[i-1]; 3555 } 3556 b->singlemalloc = PETSC_TRUE; 3557 b->free_a = PETSC_TRUE; 3558 b->free_ij = PETSC_TRUE; 3559 } else { 3560 b->free_a = PETSC_FALSE; 3561 b->free_ij = PETSC_FALSE; 3562 } 3563 3564 b->nz = 0; 3565 b->maxnz = nz; 3566 B->info.nz_unneeded = (double)b->maxnz; 3567 PetscFunctionReturn(0); 3568 } 3569 EXTERN_C_END 3570 3571 #undef __FUNCT__ 3572 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3573 /*@ 3574 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3575 3576 Input Parameters: 3577 + B - the matrix 3578 . i - the indices into j for the start of each row (starts with zero) 3579 . j - the column indices for each row (starts with zero) these must be sorted for each row 3580 - v - optional values in the matrix 3581 3582 Level: developer 3583 3584 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3585 3586 .keywords: matrix, aij, compressed row, sparse, sequential 3587 3588 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3589 @*/ 3590 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3591 { 3592 PetscErrorCode ierr; 3593 3594 PetscFunctionBegin; 3595 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3596 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3597 PetscFunctionReturn(0); 3598 } 3599 3600 EXTERN_C_BEGIN 3601 #undef __FUNCT__ 3602 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3603 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3604 { 3605 PetscInt i; 3606 PetscInt m,n; 3607 PetscInt nz; 3608 PetscInt *nnz, nz_max = 0; 3609 PetscScalar *values; 3610 PetscErrorCode ierr; 3611 3612 PetscFunctionBegin; 3613 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3614 3615 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3616 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3617 for(i = 0; i < m; i++) { 3618 nz = Ii[i+1]- Ii[i]; 3619 nz_max = PetscMax(nz_max, nz); 3620 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3621 nnz[i] = nz; 3622 } 3623 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3624 ierr = PetscFree(nnz);CHKERRQ(ierr); 3625 3626 if (v) { 3627 values = (PetscScalar*) v; 3628 } else { 3629 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3630 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3631 } 3632 3633 for(i = 0; i < m; i++) { 3634 nz = Ii[i+1] - Ii[i]; 3635 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3636 } 3637 3638 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3639 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3640 3641 if (!v) { 3642 ierr = PetscFree(values);CHKERRQ(ierr); 3643 } 3644 PetscFunctionReturn(0); 3645 } 3646 EXTERN_C_END 3647 3648 #include <../src/mat/impls/dense/seq/dense.h> 3649 #include <private/petscaxpy.h> 3650 3651 #undef __FUNCT__ 3652 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3653 /* 3654 Computes (B'*A')' since computing B*A directly is untenable 3655 3656 n p p 3657 ( ) ( ) ( ) 3658 m ( A ) * n ( B ) = m ( C ) 3659 ( ) ( ) ( ) 3660 3661 */ 3662 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3663 { 3664 PetscErrorCode ierr; 3665 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3666 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3667 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3668 PetscInt i,n,m,q,p; 3669 const PetscInt *ii,*idx; 3670 const PetscScalar *b,*a,*a_q; 3671 PetscScalar *c,*c_q; 3672 3673 PetscFunctionBegin; 3674 m = A->rmap->n; 3675 n = A->cmap->n; 3676 p = B->cmap->n; 3677 a = sub_a->v; 3678 b = sub_b->a; 3679 c = sub_c->v; 3680 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3681 3682 ii = sub_b->i; 3683 idx = sub_b->j; 3684 for (i=0; i<n; i++) { 3685 q = ii[i+1] - ii[i]; 3686 while (q-->0) { 3687 c_q = c + m*(*idx); 3688 a_q = a + m*i; 3689 PetscAXPY(c_q,*b,a_q,m); 3690 idx++; 3691 b++; 3692 } 3693 } 3694 PetscFunctionReturn(0); 3695 } 3696 3697 #undef __FUNCT__ 3698 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3699 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3700 { 3701 PetscErrorCode ierr; 3702 PetscInt m=A->rmap->n,n=B->cmap->n; 3703 Mat Cmat; 3704 3705 PetscFunctionBegin; 3706 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); 3707 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3708 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3709 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3710 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3711 Cmat->assembled = PETSC_TRUE; 3712 *C = Cmat; 3713 PetscFunctionReturn(0); 3714 } 3715 3716 /* ----------------------------------------------------------------*/ 3717 #undef __FUNCT__ 3718 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3719 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3720 { 3721 PetscErrorCode ierr; 3722 3723 PetscFunctionBegin; 3724 if (scall == MAT_INITIAL_MATRIX){ 3725 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3726 } 3727 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3728 PetscFunctionReturn(0); 3729 } 3730 3731 3732 /*MC 3733 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3734 based on compressed sparse row format. 3735 3736 Options Database Keys: 3737 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3738 3739 Level: beginner 3740 3741 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3742 M*/ 3743 3744 EXTERN_C_BEGIN 3745 #if defined(PETSC_HAVE_PASTIX) 3746 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3747 #endif 3748 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3749 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3750 #endif 3751 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3752 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3753 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3754 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3755 #if defined(PETSC_HAVE_MUMPS) 3756 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3757 #endif 3758 #if defined(PETSC_HAVE_SUPERLU) 3759 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3760 #endif 3761 #if defined(PETSC_HAVE_SUPERLU_DIST) 3762 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3763 #endif 3764 #if defined(PETSC_HAVE_SPOOLES) 3765 extern PetscErrorCode MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3766 #endif 3767 #if defined(PETSC_HAVE_UMFPACK) 3768 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3769 #endif 3770 #if defined(PETSC_HAVE_CHOLMOD) 3771 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3772 #endif 3773 #if defined(PETSC_HAVE_LUSOL) 3774 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3775 #endif 3776 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3777 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3778 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3779 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3780 #endif 3781 EXTERN_C_END 3782 3783 EXTERN_C_BEGIN 3784 #undef __FUNCT__ 3785 #define __FUNCT__ "MatCreate_SeqAIJ" 3786 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3787 { 3788 Mat_SeqAIJ *b; 3789 PetscErrorCode ierr; 3790 PetscMPIInt size; 3791 3792 PetscFunctionBegin; 3793 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3794 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3795 3796 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3797 B->data = (void*)b; 3798 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3799 b->row = 0; 3800 b->col = 0; 3801 b->icol = 0; 3802 b->reallocs = 0; 3803 b->ignorezeroentries = PETSC_FALSE; 3804 b->roworiented = PETSC_TRUE; 3805 b->nonew = 0; 3806 b->diag = 0; 3807 b->solve_work = 0; 3808 B->spptr = 0; 3809 b->saved_values = 0; 3810 b->idiag = 0; 3811 b->mdiag = 0; 3812 b->ssor_work = 0; 3813 b->omega = 1.0; 3814 b->fshift = 0.0; 3815 b->idiagvalid = PETSC_FALSE; 3816 b->ibdiagvalid = PETSC_FALSE; 3817 b->keepnonzeropattern = PETSC_FALSE; 3818 b->xtoy = 0; 3819 b->XtoY = 0; 3820 B->same_nonzero = PETSC_FALSE; 3821 3822 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3823 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3824 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3826 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3827 #endif 3828 #if defined(PETSC_HAVE_PASTIX) 3829 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3830 #endif 3831 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3832 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3833 #endif 3834 #if defined(PETSC_HAVE_SUPERLU) 3835 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3836 #endif 3837 #if defined(PETSC_HAVE_SUPERLU_DIST) 3838 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3839 #endif 3840 #if defined(PETSC_HAVE_SPOOLES) 3841 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C","MatGetFactor_seqaij_spooles",MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3842 #endif 3843 #if defined(PETSC_HAVE_MUMPS) 3844 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3845 #endif 3846 #if defined(PETSC_HAVE_UMFPACK) 3847 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3848 #endif 3849 #if defined(PETSC_HAVE_CHOLMOD) 3850 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3851 #endif 3852 #if defined(PETSC_HAVE_LUSOL) 3853 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_seqaij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3854 #endif 3855 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3856 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3857 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3858 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3859 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3860 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3861 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3862 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3863 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3864 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3865 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3866 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3867 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3868 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3869 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3870 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3871 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3872 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3873 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3874 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3875 PetscFunctionReturn(0); 3876 } 3877 EXTERN_C_END 3878 3879 #if defined(PETSC_HAVE_PTHREADCLASSES) 3880 EXTERN_C_BEGIN 3881 #undef __FUNCT__ 3882 #define __FUNCT__ "MatCreate_SeqAIJPThread" 3883 PetscErrorCode MatCreate_SeqAIJPThread(Mat B) 3884 { 3885 PetscErrorCode ierr; 3886 3887 PetscFunctionBegin; 3888 ierr = MatCreate_SeqAIJ(B); 3889 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3890 B->ops->mult = MatMult_SeqAIJPThread; 3891 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJPTHREAD);CHKERRQ(ierr); 3892 PetscFunctionReturn(0); 3893 } 3894 EXTERN_C_END 3895 #endif 3896 3897 #undef __FUNCT__ 3898 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3899 /* 3900 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3901 */ 3902 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3903 { 3904 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3905 PetscErrorCode ierr; 3906 PetscInt i,m = A->rmap->n; 3907 3908 PetscFunctionBegin; 3909 c = (Mat_SeqAIJ*)C->data; 3910 3911 C->factortype = A->factortype; 3912 c->row = 0; 3913 c->col = 0; 3914 c->icol = 0; 3915 c->reallocs = 0; 3916 3917 C->assembled = PETSC_TRUE; 3918 3919 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3920 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3921 3922 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3923 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3924 for (i=0; i<m; i++) { 3925 c->imax[i] = a->imax[i]; 3926 c->ilen[i] = a->ilen[i]; 3927 } 3928 3929 /* allocate the matrix space */ 3930 if (mallocmatspace){ 3931 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3932 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3933 c->singlemalloc = PETSC_TRUE; 3934 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3935 if (m > 0) { 3936 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3937 if (cpvalues == MAT_COPY_VALUES) { 3938 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3939 } else { 3940 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3941 } 3942 } 3943 } 3944 3945 c->ignorezeroentries = a->ignorezeroentries; 3946 c->roworiented = a->roworiented; 3947 c->nonew = a->nonew; 3948 if (a->diag) { 3949 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3950 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3951 for (i=0; i<m; i++) { 3952 c->diag[i] = a->diag[i]; 3953 } 3954 } else c->diag = 0; 3955 c->solve_work = 0; 3956 c->saved_values = 0; 3957 c->idiag = 0; 3958 c->ssor_work = 0; 3959 c->keepnonzeropattern = a->keepnonzeropattern; 3960 c->free_a = PETSC_TRUE; 3961 c->free_ij = PETSC_TRUE; 3962 c->xtoy = 0; 3963 c->XtoY = 0; 3964 3965 c->nz = a->nz; 3966 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3967 C->preallocated = PETSC_TRUE; 3968 3969 c->compressedrow.use = a->compressedrow.use; 3970 c->compressedrow.nrows = a->compressedrow.nrows; 3971 c->compressedrow.check = a->compressedrow.check; 3972 if (a->compressedrow.use){ 3973 i = a->compressedrow.nrows; 3974 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3975 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3976 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3977 } else { 3978 c->compressedrow.use = PETSC_FALSE; 3979 c->compressedrow.i = PETSC_NULL; 3980 c->compressedrow.rindex = PETSC_NULL; 3981 } 3982 C->same_nonzero = A->same_nonzero; 3983 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3984 3985 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3986 PetscFunctionReturn(0); 3987 } 3988 3989 #undef __FUNCT__ 3990 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3991 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3992 { 3993 PetscErrorCode ierr; 3994 3995 PetscFunctionBegin; 3996 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3997 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3998 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3999 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4000 PetscFunctionReturn(0); 4001 } 4002 4003 #undef __FUNCT__ 4004 #define __FUNCT__ "MatLoad_SeqAIJ" 4005 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4006 { 4007 Mat_SeqAIJ *a; 4008 PetscErrorCode ierr; 4009 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4010 int fd; 4011 PetscMPIInt size; 4012 MPI_Comm comm; 4013 PetscInt bs = 1; 4014 4015 PetscFunctionBegin; 4016 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4017 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4018 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4019 4020 ierr = PetscOptionsBegin(comm,PETSC_NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4021 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 4022 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4023 4024 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4025 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4026 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4027 M = header[1]; N = header[2]; nz = header[3]; 4028 4029 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4030 4031 /* read in row lengths */ 4032 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4033 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4034 4035 /* check if sum of rowlengths is same as nz */ 4036 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4037 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); 4038 4039 /* set global size if not set already*/ 4040 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4041 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4042 } else { 4043 /* if sizes and type are already set, check if the vector global sizes are correct */ 4044 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4045 if (rows < 0 && cols < 0){ /* user might provide local size instead of global size */ 4046 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4047 } 4048 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); 4049 } 4050 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4051 a = (Mat_SeqAIJ*)newMat->data; 4052 4053 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4054 4055 /* read in nonzero values */ 4056 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4057 4058 /* set matrix "i" values */ 4059 a->i[0] = 0; 4060 for (i=1; i<= M; i++) { 4061 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4062 a->ilen[i-1] = rowlengths[i-1]; 4063 } 4064 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4065 4066 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4067 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4068 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4069 PetscFunctionReturn(0); 4070 } 4071 4072 #undef __FUNCT__ 4073 #define __FUNCT__ "MatEqual_SeqAIJ" 4074 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4075 { 4076 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 4077 PetscErrorCode ierr; 4078 #if defined(PETSC_USE_COMPLEX) 4079 PetscInt k; 4080 #endif 4081 4082 PetscFunctionBegin; 4083 /* If the matrix dimensions are not equal,or no of nonzeros */ 4084 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4085 *flg = PETSC_FALSE; 4086 PetscFunctionReturn(0); 4087 } 4088 4089 /* if the a->i are the same */ 4090 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4091 if (!*flg) PetscFunctionReturn(0); 4092 4093 /* if a->j are the same */ 4094 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4095 if (!*flg) PetscFunctionReturn(0); 4096 4097 /* if a->a are the same */ 4098 #if defined(PETSC_USE_COMPLEX) 4099 for (k=0; k<a->nz; k++){ 4100 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 4101 *flg = PETSC_FALSE; 4102 PetscFunctionReturn(0); 4103 } 4104 } 4105 #else 4106 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4107 #endif 4108 PetscFunctionReturn(0); 4109 } 4110 4111 #undef __FUNCT__ 4112 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4113 /*@ 4114 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4115 provided by the user. 4116 4117 Collective on MPI_Comm 4118 4119 Input Parameters: 4120 + comm - must be an MPI communicator of size 1 4121 . m - number of rows 4122 . n - number of columns 4123 . i - row indices 4124 . j - column indices 4125 - a - matrix values 4126 4127 Output Parameter: 4128 . mat - the matrix 4129 4130 Level: intermediate 4131 4132 Notes: 4133 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4134 once the matrix is destroyed and not before 4135 4136 You cannot set new nonzero locations into this matrix, that will generate an error. 4137 4138 The i and j indices are 0 based 4139 4140 The format which is used for the sparse matrix input, is equivalent to a 4141 row-major ordering.. i.e for the following matrix, the input data expected is 4142 as shown: 4143 4144 1 0 0 4145 2 0 3 4146 4 5 6 4147 4148 i = {0,1,3,6} [size = nrow+1 = 3+1] 4149 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4150 v = {1,2,3,4,5,6} [size = nz = 6] 4151 4152 4153 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4154 4155 @*/ 4156 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 4157 { 4158 PetscErrorCode ierr; 4159 PetscInt ii; 4160 Mat_SeqAIJ *aij; 4161 #if defined(PETSC_USE_DEBUG) 4162 PetscInt jj; 4163 #endif 4164 4165 PetscFunctionBegin; 4166 if (i[0]) { 4167 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4168 } 4169 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4170 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4171 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4172 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4173 aij = (Mat_SeqAIJ*)(*mat)->data; 4174 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4175 4176 aij->i = i; 4177 aij->j = j; 4178 aij->a = a; 4179 aij->singlemalloc = PETSC_FALSE; 4180 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4181 aij->free_a = PETSC_FALSE; 4182 aij->free_ij = PETSC_FALSE; 4183 4184 for (ii=0; ii<m; ii++) { 4185 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4186 #if defined(PETSC_USE_DEBUG) 4187 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]); 4188 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4189 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); 4190 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); 4191 } 4192 #endif 4193 } 4194 #if defined(PETSC_USE_DEBUG) 4195 for (ii=0; ii<aij->i[m]; ii++) { 4196 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4197 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]); 4198 } 4199 #endif 4200 4201 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4202 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4203 PetscFunctionReturn(0); 4204 } 4205 4206 #undef __FUNCT__ 4207 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4208 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4209 { 4210 PetscErrorCode ierr; 4211 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4212 4213 PetscFunctionBegin; 4214 if (coloring->ctype == IS_COLORING_GLOBAL) { 4215 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4216 a->coloring = coloring; 4217 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4218 PetscInt i,*larray; 4219 ISColoring ocoloring; 4220 ISColoringValue *colors; 4221 4222 /* set coloring for diagonal portion */ 4223 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4224 for (i=0; i<A->cmap->n; i++) { 4225 larray[i] = i; 4226 } 4227 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 4228 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4229 for (i=0; i<A->cmap->n; i++) { 4230 colors[i] = coloring->colors[larray[i]]; 4231 } 4232 ierr = PetscFree(larray);CHKERRQ(ierr); 4233 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4234 a->coloring = ocoloring; 4235 } 4236 PetscFunctionReturn(0); 4237 } 4238 4239 #if defined(PETSC_HAVE_ADIC) 4240 EXTERN_C_BEGIN 4241 #include <adic/ad_utils.h> 4242 EXTERN_C_END 4243 4244 #undef __FUNCT__ 4245 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 4246 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 4247 { 4248 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4249 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 4250 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 4251 ISColoringValue *color; 4252 4253 PetscFunctionBegin; 4254 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4255 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 4256 color = a->coloring->colors; 4257 /* loop over rows */ 4258 for (i=0; i<m; i++) { 4259 nz = ii[i+1] - ii[i]; 4260 /* loop over columns putting computed value into matrix */ 4261 for (j=0; j<nz; j++) { 4262 *v++ = values[color[*jj++]]; 4263 } 4264 values += nlen; /* jump to next row of derivatives */ 4265 } 4266 PetscFunctionReturn(0); 4267 } 4268 #endif 4269 4270 #undef __FUNCT__ 4271 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4272 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4273 { 4274 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4275 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4276 MatScalar *v = a->a; 4277 PetscScalar *values = (PetscScalar *)advalues; 4278 ISColoringValue *color; 4279 4280 PetscFunctionBegin; 4281 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4282 color = a->coloring->colors; 4283 /* loop over rows */ 4284 for (i=0; i<m; i++) { 4285 nz = ii[i+1] - ii[i]; 4286 /* loop over columns putting computed value into matrix */ 4287 for (j=0; j<nz; j++) { 4288 *v++ = values[color[*jj++]]; 4289 } 4290 values += nl; /* jump to next row of derivatives */ 4291 } 4292 PetscFunctionReturn(0); 4293 } 4294 4295 /* 4296 Special version for direct calls from Fortran 4297 */ 4298 #include <private/fortranimpl.h> 4299 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4300 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4301 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4302 #define matsetvaluesseqaij_ matsetvaluesseqaij 4303 #endif 4304 4305 /* Change these macros so can be used in void function */ 4306 #undef CHKERRQ 4307 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 4308 #undef SETERRQ2 4309 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4310 4311 EXTERN_C_BEGIN 4312 #undef __FUNCT__ 4313 #define __FUNCT__ "matsetvaluesseqaij_" 4314 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4315 { 4316 Mat A = *AA; 4317 PetscInt m = *mm, n = *nn; 4318 InsertMode is = *isis; 4319 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4320 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4321 PetscInt *imax,*ai,*ailen; 4322 PetscErrorCode ierr; 4323 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4324 MatScalar *ap,value,*aa; 4325 PetscBool ignorezeroentries = a->ignorezeroentries; 4326 PetscBool roworiented = a->roworiented; 4327 4328 PetscFunctionBegin; 4329 ierr = MatPreallocated(A);CHKERRQ(ierr); 4330 imax = a->imax; 4331 ai = a->i; 4332 ailen = a->ilen; 4333 aj = a->j; 4334 aa = a->a; 4335 4336 for (k=0; k<m; k++) { /* loop over added rows */ 4337 row = im[k]; 4338 if (row < 0) continue; 4339 #if defined(PETSC_USE_DEBUG) 4340 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4341 #endif 4342 rp = aj + ai[row]; ap = aa + ai[row]; 4343 rmax = imax[row]; nrow = ailen[row]; 4344 low = 0; 4345 high = nrow; 4346 for (l=0; l<n; l++) { /* loop over added columns */ 4347 if (in[l] < 0) continue; 4348 #if defined(PETSC_USE_DEBUG) 4349 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4350 #endif 4351 col = in[l]; 4352 if (roworiented) { 4353 value = v[l + k*n]; 4354 } else { 4355 value = v[k + l*m]; 4356 } 4357 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4358 4359 if (col <= lastcol) low = 0; else high = nrow; 4360 lastcol = col; 4361 while (high-low > 5) { 4362 t = (low+high)/2; 4363 if (rp[t] > col) high = t; 4364 else low = t; 4365 } 4366 for (i=low; i<high; i++) { 4367 if (rp[i] > col) break; 4368 if (rp[i] == col) { 4369 if (is == ADD_VALUES) ap[i] += value; 4370 else ap[i] = value; 4371 goto noinsert; 4372 } 4373 } 4374 if (value == 0.0 && ignorezeroentries) goto noinsert; 4375 if (nonew == 1) goto noinsert; 4376 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4377 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4378 N = nrow++ - 1; a->nz++; high++; 4379 /* shift up all the later entries in this row */ 4380 for (ii=N; ii>=i; ii--) { 4381 rp[ii+1] = rp[ii]; 4382 ap[ii+1] = ap[ii]; 4383 } 4384 rp[i] = col; 4385 ap[i] = value; 4386 noinsert:; 4387 low = i + 1; 4388 } 4389 ailen[row] = nrow; 4390 } 4391 A->same_nonzero = PETSC_FALSE; 4392 PetscFunctionReturnVoid(); 4393 } 4394 EXTERN_C_END 4395 4396