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