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