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