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