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