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