1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the AIJ (compressed row) 5 matrix storage format. 6 */ 7 8 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 9 #include "src/inline/spops.h" 10 #include "src/inline/dot.h" 11 #include "petscbt.h" 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 15 /* only works if matrix has a full set of diagonal entries */ 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->m; 21 PetscScalar *v,*aa = aij->a; 22 23 PetscFunctionBegin; 24 ierr = MatMarkDiagonal_SeqAIJ(Y);CHKERRQ(ierr); 25 diag = aij->diag; 26 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 27 if (is == INSERT_VALUES) { 28 for (i=0; i<m; i++) { 29 aa[diag[i]] = v[i]; 30 } 31 } else { 32 for (i=0; i<m; i++) { 33 aa[diag[i]] += v[i]; 34 } 35 } 36 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 42 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 43 { 44 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 45 PetscErrorCode ierr; 46 PetscInt i,ishift; 47 48 PetscFunctionBegin; 49 *m = A->m; 50 if (!ia) PetscFunctionReturn(0); 51 ishift = 0; 52 if (symmetric && !A->structurally_symmetric) { 53 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 54 } else if (oshift == 1) { 55 PetscInt nz = a->i[A->m]; 56 /* malloc space and add 1 to i and j indices */ 57 ierr = PetscMalloc((A->m+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 58 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 59 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 60 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 61 } else { 62 *ia = a->i; *ja = a->j; 63 } 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 69 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 70 { 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 if (!ia) PetscFunctionReturn(0); 75 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 76 ierr = PetscFree(*ia);CHKERRQ(ierr); 77 ierr = PetscFree(*ja);CHKERRQ(ierr); 78 } 79 PetscFunctionReturn(0); 80 } 81 82 #undef __FUNCT__ 83 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 84 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 85 { 86 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 87 PetscErrorCode ierr; 88 PetscInt i,*collengths,*cia,*cja,n = A->n,m = A->m; 89 PetscInt nz = a->i[m],row,*jj,mr,col; 90 91 PetscFunctionBegin; 92 *nn = A->n; 93 if (!ia) PetscFunctionReturn(0); 94 if (symmetric) { 95 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 96 } else { 97 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 98 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 99 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 100 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 101 jj = a->j; 102 for (i=0; i<nz; i++) { 103 collengths[jj[i]]++; 104 } 105 cia[0] = oshift; 106 for (i=0; i<n; i++) { 107 cia[i+1] = cia[i] + collengths[i]; 108 } 109 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 110 jj = a->j; 111 for (row=0; row<m; row++) { 112 mr = a->i[row+1] - a->i[row]; 113 for (i=0; i<mr; i++) { 114 col = *jj++; 115 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 116 } 117 } 118 ierr = PetscFree(collengths);CHKERRQ(ierr); 119 *ia = cia; *ja = cja; 120 } 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 126 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 127 { 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 if (!ia) PetscFunctionReturn(0); 132 133 ierr = PetscFree(*ia);CHKERRQ(ierr); 134 ierr = PetscFree(*ja);CHKERRQ(ierr); 135 136 PetscFunctionReturn(0); 137 } 138 139 #define CHUNKSIZE 15 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "MatSetValues_SeqAIJ" 143 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 144 { 145 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 146 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 147 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 148 PetscErrorCode ierr; 149 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 150 PetscScalar *ap,value,*aa = a->a; 151 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 152 PetscTruth roworiented = a->roworiented; 153 154 PetscFunctionBegin; 155 for (k=0; k<m; k++) { /* loop over added rows */ 156 row = im[k]; 157 if (row < 0) continue; 158 #if defined(PETSC_USE_DEBUG) 159 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 160 #endif 161 rp = aj + ai[row]; ap = aa + ai[row]; 162 rmax = imax[row]; nrow = ailen[row]; 163 low = 0; 164 high = nrow; 165 for (l=0; l<n; l++) { /* loop over added columns */ 166 if (in[l] < 0) continue; 167 #if defined(PETSC_USE_DEBUG) 168 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 169 #endif 170 col = in[l]; 171 if (roworiented) { 172 value = v[l + k*n]; 173 } else { 174 value = v[k + l*m]; 175 } 176 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 177 178 if (col <= lastcol) low = 0; else high = nrow; 179 lastcol = col; 180 while (high-low > 5) { 181 t = (low+high)/2; 182 if (rp[t] > col) high = t; 183 else low = t; 184 } 185 for (i=low; i<high; i++) { 186 if (rp[i] > col) break; 187 if (rp[i] == col) { 188 if (is == ADD_VALUES) ap[i] += value; 189 else ap[i] = value; 190 goto noinsert; 191 } 192 } 193 if (value == 0.0 && ignorezeroentries) goto noinsert; 194 if (nonew == 1) goto noinsert; 195 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 196 MatSeqXAIJReallocateAIJ(a,1,nrow,row,col,rmax,aa,ai,aj,A->m,rp,ap,imax,nonew); 197 N = nrow++ - 1; a->nz++; 198 /* shift up all the later entries in this row */ 199 for (ii=N; ii>=i; ii--) { 200 rp[ii+1] = rp[ii]; 201 ap[ii+1] = ap[ii]; 202 } 203 rp[i] = col; 204 ap[i] = value; 205 noinsert:; 206 low = i + 1; 207 } 208 ailen[row] = nrow; 209 } 210 A->same_nonzero = PETSC_FALSE; 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatGetValues_SeqAIJ" 216 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 217 { 218 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 219 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 220 PetscInt *ai = a->i,*ailen = a->ilen; 221 PetscScalar *ap,*aa = a->a,zero = 0.0; 222 223 PetscFunctionBegin; 224 for (k=0; k<m; k++) { /* loop over rows */ 225 row = im[k]; 226 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 227 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 228 rp = aj + ai[row]; ap = aa + ai[row]; 229 nrow = ailen[row]; 230 for (l=0; l<n; l++) { /* loop over columns */ 231 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 232 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 233 col = in[l] ; 234 high = nrow; low = 0; /* assume unsorted */ 235 while (high-low > 5) { 236 t = (low+high)/2; 237 if (rp[t] > col) high = t; 238 else low = t; 239 } 240 for (i=low; i<high; i++) { 241 if (rp[i] > col) break; 242 if (rp[i] == col) { 243 *v++ = ap[i]; 244 goto finished; 245 } 246 } 247 *v++ = zero; 248 finished:; 249 } 250 } 251 PetscFunctionReturn(0); 252 } 253 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatView_SeqAIJ_Binary" 257 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 258 { 259 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 260 PetscErrorCode ierr; 261 PetscInt i,*col_lens; 262 int fd; 263 264 PetscFunctionBegin; 265 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 266 ierr = PetscMalloc((4+A->m)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 267 col_lens[0] = MAT_FILE_COOKIE; 268 col_lens[1] = A->m; 269 col_lens[2] = A->n; 270 col_lens[3] = a->nz; 271 272 /* store lengths of each row and write (including header) to file */ 273 for (i=0; i<A->m; i++) { 274 col_lens[4+i] = a->i[i+1] - a->i[i]; 275 } 276 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 277 ierr = PetscFree(col_lens);CHKERRQ(ierr); 278 279 /* store column indices (zero start index) */ 280 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 281 282 /* store nonzero values */ 283 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 284 PetscFunctionReturn(0); 285 } 286 287 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 291 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 292 { 293 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 294 PetscErrorCode ierr; 295 PetscInt i,j,m = A->m,shift=0; 296 const char *name; 297 PetscViewerFormat format; 298 299 PetscFunctionBegin; 300 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 301 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 302 if (format == PETSC_VIEWER_ASCII_MATLAB) { 303 PetscInt nofinalvalue = 0; 304 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 305 nofinalvalue = 1; 306 } 307 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 308 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->n);CHKERRQ(ierr); 309 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 310 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 311 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 312 313 for (i=0; i<m; i++) { 314 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 315 #if defined(PETSC_USE_COMPLEX) 316 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); 317 #else 318 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 319 #endif 320 } 321 } 322 if (nofinalvalue) { 323 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 324 } 325 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 326 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 327 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 328 PetscFunctionReturn(0); 329 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 330 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 331 for (i=0; i<m; i++) { 332 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 333 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 334 #if defined(PETSC_USE_COMPLEX) 335 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 336 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 337 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 338 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 339 } else if (PetscRealPart(a->a[j]) != 0.0) { 340 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 341 } 342 #else 343 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 344 #endif 345 } 346 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 347 } 348 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 349 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 350 PetscInt nzd=0,fshift=1,*sptr; 351 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 352 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 353 for (i=0; i<m; i++) { 354 sptr[i] = nzd+1; 355 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 356 if (a->j[j] >= i) { 357 #if defined(PETSC_USE_COMPLEX) 358 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 359 #else 360 if (a->a[j] != 0.0) nzd++; 361 #endif 362 } 363 } 364 } 365 sptr[m] = nzd+1; 366 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 367 for (i=0; i<m+1; i+=6) { 368 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);} 369 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);} 370 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);} 371 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 372 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 373 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 374 } 375 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 376 ierr = PetscFree(sptr);CHKERRQ(ierr); 377 for (i=0; i<m; i++) { 378 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 379 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 380 } 381 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 382 } 383 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 384 for (i=0; i<m; i++) { 385 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 386 if (a->j[j] >= i) { 387 #if defined(PETSC_USE_COMPLEX) 388 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 389 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 390 } 391 #else 392 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 393 #endif 394 } 395 } 396 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 397 } 398 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 399 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 400 PetscInt cnt = 0,jcnt; 401 PetscScalar value; 402 403 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 404 for (i=0; i<m; i++) { 405 jcnt = 0; 406 for (j=0; j<A->n; j++) { 407 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 408 value = a->a[cnt++]; 409 jcnt++; 410 } else { 411 value = 0.0; 412 } 413 #if defined(PETSC_USE_COMPLEX) 414 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 415 #else 416 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 417 #endif 418 } 419 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 420 } 421 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 422 } else { 423 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 424 for (i=0; i<m; i++) { 425 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 426 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 427 #if defined(PETSC_USE_COMPLEX) 428 if (PetscImaginaryPart(a->a[j]) > 0.0) { 429 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 430 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 431 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 432 } else { 433 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 434 } 435 #else 436 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 437 #endif 438 } 439 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 440 } 441 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 442 } 443 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNCT__ 448 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 449 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 450 { 451 Mat A = (Mat) Aa; 452 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 453 PetscErrorCode ierr; 454 PetscInt i,j,m = A->m,color; 455 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 456 PetscViewer viewer; 457 PetscViewerFormat format; 458 459 PetscFunctionBegin; 460 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 461 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 462 463 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 464 /* loop over matrix elements drawing boxes */ 465 466 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 467 /* Blue for negative, Cyan for zero and Red for positive */ 468 color = PETSC_DRAW_BLUE; 469 for (i=0; i<m; i++) { 470 y_l = m - i - 1.0; y_r = y_l + 1.0; 471 for (j=a->i[i]; j<a->i[i+1]; j++) { 472 x_l = a->j[j] ; x_r = x_l + 1.0; 473 #if defined(PETSC_USE_COMPLEX) 474 if (PetscRealPart(a->a[j]) >= 0.) continue; 475 #else 476 if (a->a[j] >= 0.) continue; 477 #endif 478 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 479 } 480 } 481 color = PETSC_DRAW_CYAN; 482 for (i=0; i<m; i++) { 483 y_l = m - i - 1.0; y_r = y_l + 1.0; 484 for (j=a->i[i]; j<a->i[i+1]; j++) { 485 x_l = a->j[j]; x_r = x_l + 1.0; 486 if (a->a[j] != 0.) continue; 487 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 488 } 489 } 490 color = PETSC_DRAW_RED; 491 for (i=0; i<m; i++) { 492 y_l = m - i - 1.0; y_r = y_l + 1.0; 493 for (j=a->i[i]; j<a->i[i+1]; j++) { 494 x_l = a->j[j]; x_r = x_l + 1.0; 495 #if defined(PETSC_USE_COMPLEX) 496 if (PetscRealPart(a->a[j]) <= 0.) continue; 497 #else 498 if (a->a[j] <= 0.) continue; 499 #endif 500 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 501 } 502 } 503 } else { 504 /* use contour shading to indicate magnitude of values */ 505 /* first determine max of all nonzero values */ 506 PetscInt nz = a->nz,count; 507 PetscDraw popup; 508 PetscReal scale; 509 510 for (i=0; i<nz; i++) { 511 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 512 } 513 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 514 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 515 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 516 count = 0; 517 for (i=0; i<m; i++) { 518 y_l = m - i - 1.0; y_r = y_l + 1.0; 519 for (j=a->i[i]; j<a->i[i+1]; j++) { 520 x_l = a->j[j]; x_r = x_l + 1.0; 521 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 522 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 523 count++; 524 } 525 } 526 } 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "MatView_SeqAIJ_Draw" 532 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 533 { 534 PetscErrorCode ierr; 535 PetscDraw draw; 536 PetscReal xr,yr,xl,yl,h,w; 537 PetscTruth isnull; 538 539 PetscFunctionBegin; 540 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 541 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 542 if (isnull) PetscFunctionReturn(0); 543 544 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 545 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 546 xr += w; yr += h; xl = -w; yl = -h; 547 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 548 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 549 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 553 #undef __FUNCT__ 554 #define __FUNCT__ "MatView_SeqAIJ" 555 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 556 { 557 PetscErrorCode ierr; 558 PetscTruth issocket,iascii,isbinary,isdraw; 559 560 PetscFunctionBegin; 561 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 562 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 563 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 564 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 565 if (iascii) { 566 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 567 #if defined(PETSC_USE_SOCKET_VIEWER) 568 } else if (issocket) { 569 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 570 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 571 #endif 572 } else if (isbinary) { 573 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 574 } else if (isdraw) { 575 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 576 } else { 577 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 578 } 579 ierr = MatView_Inode(A,viewer);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 585 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 586 { 587 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 588 PetscErrorCode ierr; 589 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 590 PetscInt m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 591 PetscScalar *aa = a->a,*ap; 592 PetscReal ratio=0.6; 593 594 PetscFunctionBegin; 595 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 596 597 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 598 for (i=1; i<m; i++) { 599 /* move each row back by the amount of empty slots (fshift) before it*/ 600 fshift += imax[i-1] - ailen[i-1]; 601 rmax = PetscMax(rmax,ailen[i]); 602 if (fshift) { 603 ip = aj + ai[i] ; 604 ap = aa + ai[i] ; 605 N = ailen[i]; 606 for (j=0; j<N; j++) { 607 ip[j-fshift] = ip[j]; 608 ap[j-fshift] = ap[j]; 609 } 610 } 611 ai[i] = ai[i-1] + ailen[i-1]; 612 } 613 if (m) { 614 fshift += imax[m-1] - ailen[m-1]; 615 ai[m] = ai[m-1] + ailen[m-1]; 616 } 617 /* reset ilen and imax for each row */ 618 for (i=0; i<m; i++) { 619 ailen[i] = imax[i] = ai[i+1] - ai[i]; 620 } 621 a->nz = ai[m]; 622 623 /* diagonals may have moved, so kill the diagonal pointers */ 624 if (fshift && a->diag) { 625 ierr = PetscFree(a->diag);CHKERRQ(ierr); 626 ierr = PetscLogObjectMemory(A,-(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 627 a->diag = 0; 628 } 629 ierr = PetscLogInfo((A,"MatAssemblyEnd_SeqAIJ:Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->n,fshift,a->nz));CHKERRQ(ierr); 630 ierr = PetscLogInfo((A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %D\n",a->reallocs));CHKERRQ(ierr); 631 ierr = PetscLogInfo((A,"MatAssemblyEnd_SeqAIJ:Maximum nonzeros in any row is %D\n",rmax));CHKERRQ(ierr); 632 a->reallocs = 0; 633 A->info.nz_unneeded = (double)fshift; 634 a->rmax = rmax; 635 636 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 637 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 638 A->same_nonzero = PETSC_TRUE; 639 640 ierr = MatAssemblyEnd_Inode(A,mode);CHKERRQ(ierr); 641 PetscFunctionReturn(0); 642 } 643 644 #undef __FUNCT__ 645 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 646 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 647 { 648 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 649 PetscErrorCode ierr; 650 651 PetscFunctionBegin; 652 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 #undef __FUNCT__ 657 #define __FUNCT__ "MatDestroy_SeqAIJ" 658 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 659 { 660 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 661 PetscErrorCode ierr; 662 663 PetscFunctionBegin; 664 #if defined(PETSC_USE_LOG) 665 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->m,A->n,a->nz); 666 #endif 667 if (a->freedata){ 668 ierr = MatSeqXAIJFreeAIJ(a->singlemalloc,&a->a,&a->j,&a->i);CHKERRQ(ierr); 669 } 670 if (a->row) { 671 ierr = ISDestroy(a->row);CHKERRQ(ierr); 672 } 673 if (a->col) { 674 ierr = ISDestroy(a->col);CHKERRQ(ierr); 675 } 676 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 677 if (a->ilen) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);} 678 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 679 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 680 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 681 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 682 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 683 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 684 if (a->compressedrow.use){ierr = PetscFree(a->compressedrow.i);} 685 686 ierr = MatDestroy_Inode(A);CHKERRQ(ierr); 687 688 ierr = PetscFree(a);CHKERRQ(ierr); 689 690 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 691 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 692 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 693 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 694 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 695 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 696 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 697 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 698 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 699 PetscFunctionReturn(0); 700 } 701 702 #undef __FUNCT__ 703 #define __FUNCT__ "MatCompress_SeqAIJ" 704 PetscErrorCode MatCompress_SeqAIJ(Mat A) 705 { 706 PetscFunctionBegin; 707 PetscFunctionReturn(0); 708 } 709 710 #undef __FUNCT__ 711 #define __FUNCT__ "MatSetOption_SeqAIJ" 712 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op) 713 { 714 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 715 PetscErrorCode ierr; 716 717 PetscFunctionBegin; 718 switch (op) { 719 case MAT_ROW_ORIENTED: 720 a->roworiented = PETSC_TRUE; 721 break; 722 case MAT_KEEP_ZEROED_ROWS: 723 a->keepzeroedrows = PETSC_TRUE; 724 break; 725 case MAT_COLUMN_ORIENTED: 726 a->roworiented = PETSC_FALSE; 727 break; 728 case MAT_COLUMNS_SORTED: 729 a->sorted = PETSC_TRUE; 730 break; 731 case MAT_COLUMNS_UNSORTED: 732 a->sorted = PETSC_FALSE; 733 break; 734 case MAT_NO_NEW_NONZERO_LOCATIONS: 735 a->nonew = 1; 736 break; 737 case MAT_NEW_NONZERO_LOCATION_ERR: 738 a->nonew = -1; 739 break; 740 case MAT_NEW_NONZERO_ALLOCATION_ERR: 741 a->nonew = -2; 742 break; 743 case MAT_YES_NEW_NONZERO_LOCATIONS: 744 a->nonew = 0; 745 break; 746 case MAT_IGNORE_ZERO_ENTRIES: 747 a->ignorezeroentries = PETSC_TRUE; 748 break; 749 case MAT_USE_COMPRESSEDROW: 750 a->compressedrow.use = PETSC_TRUE; 751 break; 752 case MAT_DO_NOT_USE_COMPRESSEDROW: 753 a->compressedrow.use = PETSC_FALSE; 754 break; 755 case MAT_ROWS_SORTED: 756 case MAT_ROWS_UNSORTED: 757 case MAT_YES_NEW_DIAGONALS: 758 case MAT_IGNORE_OFF_PROC_ENTRIES: 759 case MAT_USE_HASH_TABLE: 760 ierr = PetscLogInfo((A,"MatSetOption_SeqAIJ:Option ignored\n"));CHKERRQ(ierr); 761 break; 762 case MAT_NO_NEW_DIAGONALS: 763 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 764 default: 765 break; 766 } 767 ierr = MatSetOption_Inode(A,op);CHKERRQ(ierr); 768 PetscFunctionReturn(0); 769 } 770 771 #undef __FUNCT__ 772 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 773 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 774 { 775 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 776 PetscErrorCode ierr; 777 PetscInt i,j,n; 778 PetscScalar *x,zero = 0.0; 779 780 PetscFunctionBegin; 781 ierr = VecSet(v,zero);CHKERRQ(ierr); 782 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 783 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 784 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 785 for (i=0; i<A->m; i++) { 786 for (j=a->i[i]; j<a->i[i+1]; j++) { 787 if (a->j[j] == i) { 788 x[i] = a->a[j]; 789 break; 790 } 791 } 792 } 793 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 794 PetscFunctionReturn(0); 795 } 796 797 #undef __FUNCT__ 798 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 799 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 800 { 801 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 802 PetscScalar *x,*y; 803 PetscErrorCode ierr; 804 PetscInt m = A->m; 805 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 806 PetscScalar *v,alpha; 807 PetscInt n,i,*idx,*ii,*ridx=PETSC_NULL; 808 Mat_CompressedRow cprow = a->compressedrow; 809 PetscTruth usecprow = cprow.use; 810 #endif 811 812 PetscFunctionBegin; 813 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 814 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 815 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 816 817 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 818 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 819 #else 820 if (usecprow){ 821 m = cprow.nrows; 822 ii = cprow.i; 823 ridx = cprow.rindex; 824 } else { 825 ii = a->i; 826 } 827 for (i=0; i<m; i++) { 828 idx = a->j + ii[i] ; 829 v = a->a + ii[i] ; 830 n = ii[i+1] - ii[i]; 831 if (usecprow){ 832 alpha = x[ridx[i]]; 833 } else { 834 alpha = x[i]; 835 } 836 while (n-->0) {y[*idx++] += alpha * *v++;} 837 } 838 #endif 839 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 840 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 841 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 842 PetscFunctionReturn(0); 843 } 844 845 #undef __FUNCT__ 846 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 847 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 848 { 849 PetscScalar zero = 0.0; 850 PetscErrorCode ierr; 851 852 PetscFunctionBegin; 853 ierr = VecSet(yy,zero);CHKERRQ(ierr); 854 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 855 PetscFunctionReturn(0); 856 } 857 858 859 #undef __FUNCT__ 860 #define __FUNCT__ "MatMult_SeqAIJ" 861 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 862 { 863 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 864 PetscScalar *x,*y,*aa; 865 PetscErrorCode ierr; 866 PetscInt m=A->m,*aj,*ii; 867 PetscInt n,i,j,*ridx=PETSC_NULL; 868 PetscScalar sum; 869 PetscTruth usecprow=a->compressedrow.use; 870 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 871 PetscInt jrow; 872 #endif 873 874 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 875 #pragma disjoint(*x,*y,*aa) 876 #endif 877 878 PetscFunctionBegin; 879 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 880 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 881 aj = a->j; 882 aa = a->a; 883 ii = a->i; 884 if (usecprow){ /* use compressed row format */ 885 m = a->compressedrow.nrows; 886 ii = a->compressedrow.i; 887 ridx = a->compressedrow.rindex; 888 for (i=0; i<m; i++){ 889 n = ii[i+1] - ii[i]; 890 aj = a->j + ii[i]; 891 aa = a->a + ii[i]; 892 sum = 0.0; 893 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 894 y[*ridx++] = sum; 895 } 896 } else { /* do not use compressed row format */ 897 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 898 fortranmultaij_(&m,x,ii,aj,aa,y); 899 #else 900 for (i=0; i<m; i++) { 901 jrow = ii[i]; 902 n = ii[i+1] - jrow; 903 sum = 0.0; 904 for (j=0; j<n; j++) { 905 sum += aa[jrow]*x[aj[jrow]]; jrow++; 906 } 907 y[i] = sum; 908 } 909 #endif 910 } 911 ierr = PetscLogFlops(2*a->nz - m);CHKERRQ(ierr); 912 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 913 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "MatMultAdd_SeqAIJ" 919 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 920 { 921 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 922 PetscScalar *x,*y,*z,*aa; 923 PetscErrorCode ierr; 924 PetscInt m = A->m,*aj,*ii; 925 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 926 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 927 PetscScalar sum; 928 PetscTruth usecprow=a->compressedrow.use; 929 #endif 930 931 PetscFunctionBegin; 932 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 933 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 934 if (zz != yy) { 935 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 936 } else { 937 z = y; 938 } 939 940 aj = a->j; 941 aa = a->a; 942 ii = a->i; 943 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 944 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 945 #else 946 if (usecprow){ /* use compressed row format */ 947 if (zz != yy){ 948 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 949 } 950 m = a->compressedrow.nrows; 951 ii = a->compressedrow.i; 952 ridx = a->compressedrow.rindex; 953 for (i=0; i<m; i++){ 954 n = ii[i+1] - ii[i]; 955 aj = a->j + ii[i]; 956 aa = a->a + ii[i]; 957 sum = y[*ridx]; 958 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 959 z[*ridx++] = sum; 960 } 961 } else { /* do not use compressed row format */ 962 for (i=0; i<m; i++) { 963 jrow = ii[i]; 964 n = ii[i+1] - jrow; 965 sum = y[i]; 966 for (j=0; j<n; j++) { 967 sum += aa[jrow]*x[aj[jrow]]; jrow++; 968 } 969 z[i] = sum; 970 } 971 } 972 #endif 973 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 974 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 975 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 976 if (zz != yy) { 977 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 978 } 979 PetscFunctionReturn(0); 980 } 981 982 /* 983 Adds diagonal pointers to sparse matrix structure. 984 */ 985 #undef __FUNCT__ 986 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 987 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 988 { 989 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 990 PetscErrorCode ierr; 991 PetscInt i,j,*diag,m = A->m; 992 993 PetscFunctionBegin; 994 if (a->diag) PetscFunctionReturn(0); 995 996 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&diag);CHKERRQ(ierr); 997 ierr = PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 998 for (i=0; i<A->m; i++) { 999 diag[i] = a->i[i+1]; 1000 for (j=a->i[i]; j<a->i[i+1]; j++) { 1001 if (a->j[j] == i) { 1002 diag[i] = j; 1003 break; 1004 } 1005 } 1006 } 1007 a->diag = diag; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 /* 1012 Checks for missing diagonals 1013 */ 1014 #undef __FUNCT__ 1015 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1016 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A) 1017 { 1018 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1019 PetscErrorCode ierr; 1020 PetscInt *diag,*jj = a->j,i; 1021 1022 PetscFunctionBegin; 1023 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1024 diag = a->diag; 1025 for (i=0; i<A->m; i++) { 1026 if (jj[diag[i]] != i) { 1027 SETERRQ1(PETSC_ERR_PLIB,"Matrix is missing diagonal number %D",i); 1028 } 1029 } 1030 PetscFunctionReturn(0); 1031 } 1032 1033 #undef __FUNCT__ 1034 #define __FUNCT__ "MatRelax_SeqAIJ" 1035 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1036 { 1037 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1038 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1039 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1040 PetscErrorCode ierr; 1041 PetscInt n = A->n,m = A->m,i; 1042 const PetscInt *idx,*diag; 1043 1044 PetscFunctionBegin; 1045 its = its*lits; 1046 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1047 1048 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1049 diag = a->diag; 1050 if (!a->idiag) { 1051 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1052 a->ssor = a->idiag + m; 1053 mdiag = a->ssor + m; 1054 1055 v = a->a; 1056 1057 /* this is wrong when fshift omega changes each iteration */ 1058 if (omega == 1.0 && !fshift) { 1059 for (i=0; i<m; i++) { 1060 mdiag[i] = v[diag[i]]; 1061 a->idiag[i] = 1.0/v[diag[i]]; 1062 } 1063 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1064 } else { 1065 for (i=0; i<m; i++) { 1066 mdiag[i] = v[diag[i]]; 1067 a->idiag[i] = omega/(fshift + v[diag[i]]); 1068 } 1069 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1070 } 1071 } 1072 t = a->ssor; 1073 idiag = a->idiag; 1074 mdiag = a->idiag + 2*m; 1075 1076 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1077 if (xx != bb) { 1078 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1079 } else { 1080 b = x; 1081 } 1082 1083 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1084 xs = x; 1085 if (flag == SOR_APPLY_UPPER) { 1086 /* apply (U + D/omega) to the vector */ 1087 bs = b; 1088 for (i=0; i<m; i++) { 1089 d = fshift + a->a[diag[i]]; 1090 n = a->i[i+1] - diag[i] - 1; 1091 idx = a->j + diag[i] + 1; 1092 v = a->a + diag[i] + 1; 1093 sum = b[i]*d/omega; 1094 SPARSEDENSEDOT(sum,bs,v,idx,n); 1095 x[i] = sum; 1096 } 1097 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1098 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1099 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1100 PetscFunctionReturn(0); 1101 } 1102 1103 1104 /* Let A = L + U + D; where L is lower trianglar, 1105 U is upper triangular, E is diagonal; This routine applies 1106 1107 (L + E)^{-1} A (U + E)^{-1} 1108 1109 to a vector efficiently using Eisenstat's trick. This is for 1110 the case of SSOR preconditioner, so E is D/omega where omega 1111 is the relaxation factor. 1112 */ 1113 1114 if (flag == SOR_APPLY_LOWER) { 1115 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1116 } else if (flag & SOR_EISENSTAT) { 1117 /* Let A = L + U + D; where L is lower trianglar, 1118 U is upper triangular, E is diagonal; This routine applies 1119 1120 (L + E)^{-1} A (U + E)^{-1} 1121 1122 to a vector efficiently using Eisenstat's trick. This is for 1123 the case of SSOR preconditioner, so E is D/omega where omega 1124 is the relaxation factor. 1125 */ 1126 scale = (2.0/omega) - 1.0; 1127 1128 /* x = (E + U)^{-1} b */ 1129 for (i=m-1; i>=0; i--) { 1130 n = a->i[i+1] - diag[i] - 1; 1131 idx = a->j + diag[i] + 1; 1132 v = a->a + diag[i] + 1; 1133 sum = b[i]; 1134 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1135 x[i] = sum*idiag[i]; 1136 } 1137 1138 /* t = b - (2*E - D)x */ 1139 v = a->a; 1140 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1141 1142 /* t = (E + L)^{-1}t */ 1143 ts = t; 1144 diag = a->diag; 1145 for (i=0; i<m; i++) { 1146 n = diag[i] - a->i[i]; 1147 idx = a->j + a->i[i]; 1148 v = a->a + a->i[i]; 1149 sum = t[i]; 1150 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1151 t[i] = sum*idiag[i]; 1152 /* x = x + t */ 1153 x[i] += t[i]; 1154 } 1155 1156 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1157 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1158 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1159 PetscFunctionReturn(0); 1160 } 1161 if (flag & SOR_ZERO_INITIAL_GUESS) { 1162 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1163 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1164 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1165 #else 1166 for (i=0; i<m; i++) { 1167 n = diag[i] - a->i[i]; 1168 idx = a->j + a->i[i]; 1169 v = a->a + a->i[i]; 1170 sum = b[i]; 1171 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1172 x[i] = sum*idiag[i]; 1173 } 1174 #endif 1175 xb = x; 1176 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1177 } else xb = b; 1178 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1179 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1180 for (i=0; i<m; i++) { 1181 x[i] *= mdiag[i]; 1182 } 1183 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1184 } 1185 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1186 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1187 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1188 #else 1189 for (i=m-1; i>=0; i--) { 1190 n = a->i[i+1] - diag[i] - 1; 1191 idx = a->j + diag[i] + 1; 1192 v = a->a + diag[i] + 1; 1193 sum = xb[i]; 1194 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1195 x[i] = sum*idiag[i]; 1196 } 1197 #endif 1198 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1199 } 1200 its--; 1201 } 1202 while (its--) { 1203 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1204 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1205 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1206 #else 1207 for (i=0; i<m; i++) { 1208 d = fshift + a->a[diag[i]]; 1209 n = a->i[i+1] - a->i[i]; 1210 idx = a->j + a->i[i]; 1211 v = a->a + a->i[i]; 1212 sum = b[i]; 1213 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1214 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1215 } 1216 #endif 1217 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1218 } 1219 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1220 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1221 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1222 #else 1223 for (i=m-1; i>=0; i--) { 1224 d = fshift + a->a[diag[i]]; 1225 n = a->i[i+1] - a->i[i]; 1226 idx = a->j + a->i[i]; 1227 v = a->a + a->i[i]; 1228 sum = b[i]; 1229 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1230 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1231 } 1232 #endif 1233 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1234 } 1235 } 1236 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1237 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1238 PetscFunctionReturn(0); 1239 } 1240 1241 #undef __FUNCT__ 1242 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1243 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1244 { 1245 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1246 1247 PetscFunctionBegin; 1248 info->rows_global = (double)A->m; 1249 info->columns_global = (double)A->n; 1250 info->rows_local = (double)A->m; 1251 info->columns_local = (double)A->n; 1252 info->block_size = 1.0; 1253 info->nz_allocated = (double)a->maxnz; 1254 info->nz_used = (double)a->nz; 1255 info->nz_unneeded = (double)(a->maxnz - a->nz); 1256 info->assemblies = (double)A->num_ass; 1257 info->mallocs = (double)a->reallocs; 1258 info->memory = A->mem; 1259 if (A->factor) { 1260 info->fill_ratio_given = A->info.fill_ratio_given; 1261 info->fill_ratio_needed = A->info.fill_ratio_needed; 1262 info->factor_mallocs = A->info.factor_mallocs; 1263 } else { 1264 info->fill_ratio_given = 0; 1265 info->fill_ratio_needed = 0; 1266 info->factor_mallocs = 0; 1267 } 1268 PetscFunctionReturn(0); 1269 } 1270 1271 #undef __FUNCT__ 1272 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1273 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1274 { 1275 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1276 PetscInt i,m = A->m - 1; 1277 PetscErrorCode ierr; 1278 1279 PetscFunctionBegin; 1280 if (a->keepzeroedrows) { 1281 for (i=0; i<N; i++) { 1282 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1283 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1284 } 1285 if (diag != 0.0) { 1286 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1287 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1288 for (i=0; i<N; i++) { 1289 a->a[a->diag[rows[i]]] = diag; 1290 } 1291 } 1292 A->same_nonzero = PETSC_TRUE; 1293 } else { 1294 if (diag != 0.0) { 1295 for (i=0; i<N; i++) { 1296 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1297 if (a->ilen[rows[i]] > 0) { 1298 a->ilen[rows[i]] = 1; 1299 a->a[a->i[rows[i]]] = diag; 1300 a->j[a->i[rows[i]]] = rows[i]; 1301 } else { /* in case row was completely empty */ 1302 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1303 } 1304 } 1305 } else { 1306 for (i=0; i<N; i++) { 1307 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1308 a->ilen[rows[i]] = 0; 1309 } 1310 } 1311 A->same_nonzero = PETSC_FALSE; 1312 } 1313 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1314 PetscFunctionReturn(0); 1315 } 1316 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "MatGetRow_SeqAIJ" 1319 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1320 { 1321 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1322 PetscInt *itmp; 1323 1324 PetscFunctionBegin; 1325 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1326 1327 *nz = a->i[row+1] - a->i[row]; 1328 if (v) *v = a->a + a->i[row]; 1329 if (idx) { 1330 itmp = a->j + a->i[row]; 1331 if (*nz) { 1332 *idx = itmp; 1333 } 1334 else *idx = 0; 1335 } 1336 PetscFunctionReturn(0); 1337 } 1338 1339 /* remove this function? */ 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1342 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1343 { 1344 PetscFunctionBegin; 1345 PetscFunctionReturn(0); 1346 } 1347 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "MatNorm_SeqAIJ" 1350 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1351 { 1352 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1353 PetscScalar *v = a->a; 1354 PetscReal sum = 0.0; 1355 PetscErrorCode ierr; 1356 PetscInt i,j; 1357 1358 PetscFunctionBegin; 1359 if (type == NORM_FROBENIUS) { 1360 for (i=0; i<a->nz; i++) { 1361 #if defined(PETSC_USE_COMPLEX) 1362 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1363 #else 1364 sum += (*v)*(*v); v++; 1365 #endif 1366 } 1367 *nrm = sqrt(sum); 1368 } else if (type == NORM_1) { 1369 PetscReal *tmp; 1370 PetscInt *jj = a->j; 1371 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1372 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1373 *nrm = 0.0; 1374 for (j=0; j<a->nz; j++) { 1375 tmp[*jj++] += PetscAbsScalar(*v); v++; 1376 } 1377 for (j=0; j<A->n; j++) { 1378 if (tmp[j] > *nrm) *nrm = tmp[j]; 1379 } 1380 ierr = PetscFree(tmp);CHKERRQ(ierr); 1381 } else if (type == NORM_INFINITY) { 1382 *nrm = 0.0; 1383 for (j=0; j<A->m; j++) { 1384 v = a->a + a->i[j]; 1385 sum = 0.0; 1386 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1387 sum += PetscAbsScalar(*v); v++; 1388 } 1389 if (sum > *nrm) *nrm = sum; 1390 } 1391 } else { 1392 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1393 } 1394 PetscFunctionReturn(0); 1395 } 1396 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "MatTranspose_SeqAIJ" 1399 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1400 { 1401 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1402 Mat C; 1403 PetscErrorCode ierr; 1404 PetscInt i,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1405 PetscScalar *array = a->a; 1406 1407 PetscFunctionBegin; 1408 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1409 ierr = PetscMalloc((1+A->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1410 ierr = PetscMemzero(col,(1+A->n)*sizeof(PetscInt));CHKERRQ(ierr); 1411 1412 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1413 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1414 ierr = MatSetSizes(C,A->n,m,A->n,m);CHKERRQ(ierr); 1415 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1416 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1417 ierr = PetscFree(col);CHKERRQ(ierr); 1418 for (i=0; i<m; i++) { 1419 len = ai[i+1]-ai[i]; 1420 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1421 array += len; 1422 aj += len; 1423 } 1424 1425 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1426 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1427 1428 if (B) { 1429 *B = C; 1430 } else { 1431 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1432 } 1433 PetscFunctionReturn(0); 1434 } 1435 1436 EXTERN_C_BEGIN 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1439 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1440 { 1441 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1442 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1443 PetscErrorCode ierr; 1444 PetscInt ma,na,mb,nb, i; 1445 1446 PetscFunctionBegin; 1447 bij = (Mat_SeqAIJ *) B->data; 1448 1449 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1450 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1451 if (ma!=nb || na!=mb){ 1452 *f = PETSC_FALSE; 1453 PetscFunctionReturn(0); 1454 } 1455 aii = aij->i; bii = bij->i; 1456 adx = aij->j; bdx = bij->j; 1457 va = aij->a; vb = bij->a; 1458 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1459 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1460 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1461 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1462 1463 *f = PETSC_TRUE; 1464 for (i=0; i<ma; i++) { 1465 while (aptr[i]<aii[i+1]) { 1466 PetscInt idc,idr; 1467 PetscScalar vc,vr; 1468 /* column/row index/value */ 1469 idc = adx[aptr[i]]; 1470 idr = bdx[bptr[idc]]; 1471 vc = va[aptr[i]]; 1472 vr = vb[bptr[idc]]; 1473 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1474 *f = PETSC_FALSE; 1475 goto done; 1476 } else { 1477 aptr[i]++; 1478 if (B || i!=idc) bptr[idc]++; 1479 } 1480 } 1481 } 1482 done: 1483 ierr = PetscFree(aptr);CHKERRQ(ierr); 1484 if (B) { 1485 ierr = PetscFree(bptr);CHKERRQ(ierr); 1486 } 1487 PetscFunctionReturn(0); 1488 } 1489 EXTERN_C_END 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1493 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1494 { 1495 PetscErrorCode ierr; 1496 PetscFunctionBegin; 1497 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1498 PetscFunctionReturn(0); 1499 } 1500 1501 #undef __FUNCT__ 1502 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1503 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1504 { 1505 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1506 PetscScalar *l,*r,x,*v; 1507 PetscErrorCode ierr; 1508 PetscInt i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1509 1510 PetscFunctionBegin; 1511 if (ll) { 1512 /* The local size is used so that VecMPI can be passed to this routine 1513 by MatDiagonalScale_MPIAIJ */ 1514 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1515 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1516 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1517 v = a->a; 1518 for (i=0; i<m; i++) { 1519 x = l[i]; 1520 M = a->i[i+1] - a->i[i]; 1521 for (j=0; j<M; j++) { (*v++) *= x;} 1522 } 1523 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1524 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1525 } 1526 if (rr) { 1527 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1528 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1529 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1530 v = a->a; jj = a->j; 1531 for (i=0; i<nz; i++) { 1532 (*v++) *= r[*jj++]; 1533 } 1534 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1535 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1536 } 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1542 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1543 { 1544 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1545 PetscErrorCode ierr; 1546 PetscInt *smap,i,k,kstart,kend,oldcols = A->n,*lens; 1547 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1548 PetscInt *irow,*icol,nrows,ncols; 1549 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1550 PetscScalar *a_new,*mat_a; 1551 Mat C; 1552 PetscTruth stride; 1553 1554 PetscFunctionBegin; 1555 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1556 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1557 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1558 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1559 1560 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1561 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1562 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1563 1564 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1565 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1566 if (stride && step == 1) { 1567 /* special case of contiguous rows */ 1568 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1569 starts = lens + nrows; 1570 /* loop over new rows determining lens and starting points */ 1571 for (i=0; i<nrows; i++) { 1572 kstart = ai[irow[i]]; 1573 kend = kstart + ailen[irow[i]]; 1574 for (k=kstart; k<kend; k++) { 1575 if (aj[k] >= first) { 1576 starts[i] = k; 1577 break; 1578 } 1579 } 1580 sum = 0; 1581 while (k < kend) { 1582 if (aj[k++] >= first+ncols) break; 1583 sum++; 1584 } 1585 lens[i] = sum; 1586 } 1587 /* create submatrix */ 1588 if (scall == MAT_REUSE_MATRIX) { 1589 PetscInt n_cols,n_rows; 1590 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1591 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1592 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1593 C = *B; 1594 } else { 1595 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1596 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1597 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1598 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1599 } 1600 c = (Mat_SeqAIJ*)C->data; 1601 1602 /* loop over rows inserting into submatrix */ 1603 a_new = c->a; 1604 j_new = c->j; 1605 i_new = c->i; 1606 1607 for (i=0; i<nrows; i++) { 1608 ii = starts[i]; 1609 lensi = lens[i]; 1610 for (k=0; k<lensi; k++) { 1611 *j_new++ = aj[ii+k] - first; 1612 } 1613 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1614 a_new += lensi; 1615 i_new[i+1] = i_new[i] + lensi; 1616 c->ilen[i] = lensi; 1617 } 1618 ierr = PetscFree(lens);CHKERRQ(ierr); 1619 } else { 1620 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1621 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1622 1623 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1624 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1625 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1626 /* determine lens of each row */ 1627 for (i=0; i<nrows; i++) { 1628 kstart = ai[irow[i]]; 1629 kend = kstart + a->ilen[irow[i]]; 1630 lens[i] = 0; 1631 for (k=kstart; k<kend; k++) { 1632 if (smap[aj[k]]) { 1633 lens[i]++; 1634 } 1635 } 1636 } 1637 /* Create and fill new matrix */ 1638 if (scall == MAT_REUSE_MATRIX) { 1639 PetscTruth equal; 1640 1641 c = (Mat_SeqAIJ *)((*B)->data); 1642 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1643 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1644 if (!equal) { 1645 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1646 } 1647 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(PetscInt));CHKERRQ(ierr); 1648 C = *B; 1649 } else { 1650 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1651 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1652 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1653 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1654 } 1655 c = (Mat_SeqAIJ *)(C->data); 1656 for (i=0; i<nrows; i++) { 1657 row = irow[i]; 1658 kstart = ai[row]; 1659 kend = kstart + a->ilen[row]; 1660 mat_i = c->i[i]; 1661 mat_j = c->j + mat_i; 1662 mat_a = c->a + mat_i; 1663 mat_ilen = c->ilen + i; 1664 for (k=kstart; k<kend; k++) { 1665 if ((tcol=smap[a->j[k]])) { 1666 *mat_j++ = tcol - 1; 1667 *mat_a++ = a->a[k]; 1668 (*mat_ilen)++; 1669 1670 } 1671 } 1672 } 1673 /* Free work space */ 1674 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1675 ierr = PetscFree(smap);CHKERRQ(ierr); 1676 ierr = PetscFree(lens);CHKERRQ(ierr); 1677 } 1678 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1679 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1680 1681 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1682 *B = C; 1683 PetscFunctionReturn(0); 1684 } 1685 1686 /* 1687 */ 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1690 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1691 { 1692 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1693 PetscErrorCode ierr; 1694 Mat outA; 1695 PetscTruth row_identity,col_identity; 1696 1697 PetscFunctionBegin; 1698 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1699 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1700 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1701 if (!row_identity || !col_identity) { 1702 SETERRQ(PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU"); 1703 } 1704 1705 outA = inA; 1706 inA->factor = FACTOR_LU; 1707 a->row = row; 1708 a->col = col; 1709 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1710 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1711 1712 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1713 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1714 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1715 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1716 1717 if (!a->solve_work) { /* this matrix may have been factored before */ 1718 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1719 } 1720 1721 if (!a->diag) { 1722 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1723 } 1724 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 #include "petscblaslapack.h" 1729 #undef __FUNCT__ 1730 #define __FUNCT__ "MatScale_SeqAIJ" 1731 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1732 { 1733 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1734 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1735 PetscScalar oalpha = alpha; 1736 PetscErrorCode ierr; 1737 1738 1739 PetscFunctionBegin; 1740 BLASscal_(&bnz,&oalpha,a->a,&one); 1741 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 #undef __FUNCT__ 1746 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1747 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1748 { 1749 PetscErrorCode ierr; 1750 PetscInt i; 1751 1752 PetscFunctionBegin; 1753 if (scall == MAT_INITIAL_MATRIX) { 1754 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1755 } 1756 1757 for (i=0; i<n; i++) { 1758 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1759 } 1760 PetscFunctionReturn(0); 1761 } 1762 1763 #undef __FUNCT__ 1764 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1765 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1766 { 1767 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1768 PetscErrorCode ierr; 1769 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1770 PetscInt start,end,*ai,*aj; 1771 PetscBT table; 1772 1773 PetscFunctionBegin; 1774 m = A->m; 1775 ai = a->i; 1776 aj = a->j; 1777 1778 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1779 1780 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1781 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1782 1783 for (i=0; i<is_max; i++) { 1784 /* Initialize the two local arrays */ 1785 isz = 0; 1786 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1787 1788 /* Extract the indices, assume there can be duplicate entries */ 1789 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1790 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1791 1792 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1793 for (j=0; j<n ; ++j){ 1794 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1795 } 1796 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1797 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1798 1799 k = 0; 1800 for (j=0; j<ov; j++){ /* for each overlap */ 1801 n = isz; 1802 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1803 row = nidx[k]; 1804 start = ai[row]; 1805 end = ai[row+1]; 1806 for (l = start; l<end ; l++){ 1807 val = aj[l] ; 1808 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1809 } 1810 } 1811 } 1812 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1813 } 1814 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1815 ierr = PetscFree(nidx);CHKERRQ(ierr); 1816 PetscFunctionReturn(0); 1817 } 1818 1819 /* -------------------------------------------------------------- */ 1820 #undef __FUNCT__ 1821 #define __FUNCT__ "MatPermute_SeqAIJ" 1822 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1823 { 1824 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1825 PetscErrorCode ierr; 1826 PetscInt i,nz,m = A->m,n = A->n,*col; 1827 PetscInt *row,*cnew,j,*lens; 1828 IS icolp,irowp; 1829 PetscInt *cwork; 1830 PetscScalar *vwork; 1831 1832 PetscFunctionBegin; 1833 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1834 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1835 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1836 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1837 1838 /* determine lengths of permuted rows */ 1839 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1840 for (i=0; i<m; i++) { 1841 lens[row[i]] = a->i[i+1] - a->i[i]; 1842 } 1843 ierr = MatCreate(A->comm,B);CHKERRQ(ierr); 1844 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1845 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1846 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1847 ierr = PetscFree(lens);CHKERRQ(ierr); 1848 1849 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1850 for (i=0; i<m; i++) { 1851 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1852 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1853 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1854 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1855 } 1856 ierr = PetscFree(cnew);CHKERRQ(ierr); 1857 (*B)->assembled = PETSC_FALSE; 1858 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1859 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1860 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1861 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1862 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1863 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1864 PetscFunctionReturn(0); 1865 } 1866 1867 #undef __FUNCT__ 1868 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1869 PetscErrorCode MatPrintHelp_SeqAIJ(Mat A) 1870 { 1871 static PetscTruth called = PETSC_FALSE; 1872 MPI_Comm comm = A->comm; 1873 PetscErrorCode ierr; 1874 1875 PetscFunctionBegin; 1876 ierr = MatPrintHelp_Inode(A);CHKERRQ(ierr); 1877 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1878 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1879 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1880 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1881 ierr = (*PetscHelpPrintf)(comm," -mat_no_compressedrow: Do not use compressedrow\n");CHKERRQ(ierr); 1882 PetscFunctionReturn(0); 1883 } 1884 1885 #undef __FUNCT__ 1886 #define __FUNCT__ "MatCopy_SeqAIJ" 1887 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1888 { 1889 PetscErrorCode ierr; 1890 1891 PetscFunctionBegin; 1892 /* If the two matrices have the same copy implementation, use fast copy. */ 1893 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1894 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1895 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1896 1897 if (a->i[A->m] != b->i[B->m]) { 1898 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1899 } 1900 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1901 } else { 1902 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1903 } 1904 PetscFunctionReturn(0); 1905 } 1906 1907 #undef __FUNCT__ 1908 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1909 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1910 { 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1915 PetscFunctionReturn(0); 1916 } 1917 1918 #undef __FUNCT__ 1919 #define __FUNCT__ "MatGetArray_SeqAIJ" 1920 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1921 { 1922 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1923 PetscFunctionBegin; 1924 *array = a->a; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1930 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1931 { 1932 PetscFunctionBegin; 1933 PetscFunctionReturn(0); 1934 } 1935 1936 #undef __FUNCT__ 1937 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1938 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1939 { 1940 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 1941 PetscErrorCode ierr; 1942 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1943 PetscScalar dx,*y,*xx,*w3_array; 1944 PetscScalar *vscale_array; 1945 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1946 Vec w1,w2,w3; 1947 void *fctx = coloring->fctx; 1948 PetscTruth flg; 1949 1950 PetscFunctionBegin; 1951 if (!coloring->w1) { 1952 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1953 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 1954 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1955 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 1956 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1957 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 1958 } 1959 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1960 1961 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1962 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1963 if (flg) { 1964 ierr = PetscLogInfo((coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"));CHKERRQ(ierr); 1965 } else { 1966 PetscTruth assembled; 1967 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 1968 if (assembled) { 1969 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1970 } 1971 } 1972 1973 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1974 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1975 1976 /* 1977 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1978 coloring->F for the coarser grids from the finest 1979 */ 1980 if (coloring->F) { 1981 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1982 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1983 if (m1 != m2) { 1984 coloring->F = 0; 1985 } 1986 } 1987 1988 if (coloring->F) { 1989 w1 = coloring->F; 1990 coloring->F = 0; 1991 } else { 1992 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1993 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1994 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1995 } 1996 1997 /* 1998 Compute all the scale factors and share with other processors 1999 */ 2000 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2001 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2002 for (k=0; k<coloring->ncolors; k++) { 2003 /* 2004 Loop over each column associated with color adding the 2005 perturbation to the vector w3. 2006 */ 2007 for (l=0; l<coloring->ncolumns[k]; l++) { 2008 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2009 dx = xx[col]; 2010 if (dx == 0.0) dx = 1.0; 2011 #if !defined(PETSC_USE_COMPLEX) 2012 if (dx < umin && dx >= 0.0) dx = umin; 2013 else if (dx < 0.0 && dx > -umin) dx = -umin; 2014 #else 2015 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2016 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2017 #endif 2018 dx *= epsilon; 2019 vscale_array[col] = 1.0/dx; 2020 } 2021 } 2022 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2023 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2024 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2025 2026 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2027 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2028 2029 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2030 else vscaleforrow = coloring->columnsforrow; 2031 2032 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2033 /* 2034 Loop over each color 2035 */ 2036 for (k=0; k<coloring->ncolors; k++) { 2037 coloring->currentcolor = k; 2038 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2039 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2040 /* 2041 Loop over each column associated with color adding the 2042 perturbation to the vector w3. 2043 */ 2044 for (l=0; l<coloring->ncolumns[k]; l++) { 2045 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2046 dx = xx[col]; 2047 if (dx == 0.0) dx = 1.0; 2048 #if !defined(PETSC_USE_COMPLEX) 2049 if (dx < umin && dx >= 0.0) dx = umin; 2050 else if (dx < 0.0 && dx > -umin) dx = -umin; 2051 #else 2052 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2053 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2054 #endif 2055 dx *= epsilon; 2056 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2057 w3_array[col] += dx; 2058 } 2059 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2060 2061 /* 2062 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2063 */ 2064 2065 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2066 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2067 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2068 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2069 2070 /* 2071 Loop over rows of vector, putting results into Jacobian matrix 2072 */ 2073 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2074 for (l=0; l<coloring->nrows[k]; l++) { 2075 row = coloring->rows[k][l]; 2076 col = coloring->columnsforrow[k][l]; 2077 y[row] *= vscale_array[vscaleforrow[k][l]]; 2078 srow = row + start; 2079 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2080 } 2081 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2082 } 2083 coloring->currentcolor = k; 2084 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2085 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2086 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2087 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2088 PetscFunctionReturn(0); 2089 } 2090 2091 #include "petscblaslapack.h" 2092 #undef __FUNCT__ 2093 #define __FUNCT__ "MatAXPY_SeqAIJ" 2094 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2095 { 2096 PetscErrorCode ierr; 2097 PetscInt i; 2098 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2099 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2100 2101 PetscFunctionBegin; 2102 if (str == SAME_NONZERO_PATTERN) { 2103 PetscScalar alpha = a; 2104 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2105 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2106 if (y->xtoy && y->XtoY != X) { 2107 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2108 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2109 } 2110 if (!y->xtoy) { /* get xtoy */ 2111 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2112 y->XtoY = X; 2113 } 2114 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2115 ierr = PetscLogInfo((0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz));CHKERRQ(ierr); 2116 } else { 2117 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2118 } 2119 PetscFunctionReturn(0); 2120 } 2121 2122 #undef __FUNCT__ 2123 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2124 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2125 { 2126 PetscFunctionBegin; 2127 PetscFunctionReturn(0); 2128 } 2129 2130 #undef __FUNCT__ 2131 #define __FUNCT__ "MatConjugate_SeqAIJ" 2132 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2133 { 2134 #if defined(PETSC_USE_COMPLEX) 2135 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2136 PetscInt i,nz; 2137 PetscScalar *a; 2138 2139 PetscFunctionBegin; 2140 nz = aij->nz; 2141 a = aij->a; 2142 for (i=0; i<nz; i++) { 2143 a[i] = PetscConj(a[i]); 2144 } 2145 #else 2146 PetscFunctionBegin; 2147 #endif 2148 PetscFunctionReturn(0); 2149 } 2150 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2153 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v) 2154 { 2155 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2156 PetscErrorCode ierr; 2157 PetscInt i,j,m = A->m,*ai,*aj,ncols,n; 2158 PetscReal atmp; 2159 PetscScalar *x,zero = 0.0; 2160 MatScalar *aa; 2161 2162 PetscFunctionBegin; 2163 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2164 aa = a->a; 2165 ai = a->i; 2166 aj = a->j; 2167 2168 ierr = VecSet(v,zero);CHKERRQ(ierr); 2169 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2170 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2171 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2172 for (i=0; i<m; i++) { 2173 ncols = ai[1] - ai[0]; ai++; 2174 for (j=0; j<ncols; j++){ 2175 atmp = PetscAbsScalar(*aa); aa++; 2176 if (PetscAbsScalar(x[i]) < atmp) x[i] = atmp; 2177 aj++; 2178 } 2179 } 2180 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 /* -------------------------------------------------------------------*/ 2185 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2186 MatGetRow_SeqAIJ, 2187 MatRestoreRow_SeqAIJ, 2188 MatMult_SeqAIJ, 2189 /* 4*/ MatMultAdd_SeqAIJ, 2190 MatMultTranspose_SeqAIJ, 2191 MatMultTransposeAdd_SeqAIJ, 2192 MatSolve_SeqAIJ, 2193 MatSolveAdd_SeqAIJ, 2194 MatSolveTranspose_SeqAIJ, 2195 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2196 MatLUFactor_SeqAIJ, 2197 0, 2198 MatRelax_SeqAIJ, 2199 MatTranspose_SeqAIJ, 2200 /*15*/ MatGetInfo_SeqAIJ, 2201 MatEqual_SeqAIJ, 2202 MatGetDiagonal_SeqAIJ, 2203 MatDiagonalScale_SeqAIJ, 2204 MatNorm_SeqAIJ, 2205 /*20*/ 0, 2206 MatAssemblyEnd_SeqAIJ, 2207 MatCompress_SeqAIJ, 2208 MatSetOption_SeqAIJ, 2209 MatZeroEntries_SeqAIJ, 2210 /*25*/ MatZeroRows_SeqAIJ, 2211 MatLUFactorSymbolic_SeqAIJ, 2212 MatLUFactorNumeric_SeqAIJ, 2213 MatCholeskyFactorSymbolic_SeqAIJ, 2214 MatCholeskyFactorNumeric_SeqAIJ, 2215 /*30*/ MatSetUpPreallocation_SeqAIJ, 2216 MatILUFactorSymbolic_SeqAIJ, 2217 MatICCFactorSymbolic_SeqAIJ, 2218 MatGetArray_SeqAIJ, 2219 MatRestoreArray_SeqAIJ, 2220 /*35*/ MatDuplicate_SeqAIJ, 2221 0, 2222 0, 2223 MatILUFactor_SeqAIJ, 2224 0, 2225 /*40*/ MatAXPY_SeqAIJ, 2226 MatGetSubMatrices_SeqAIJ, 2227 MatIncreaseOverlap_SeqAIJ, 2228 MatGetValues_SeqAIJ, 2229 MatCopy_SeqAIJ, 2230 /*45*/ MatPrintHelp_SeqAIJ, 2231 MatScale_SeqAIJ, 2232 0, 2233 MatDiagonalSet_SeqAIJ, 2234 MatILUDTFactor_SeqAIJ, 2235 /*50*/ MatSetBlockSize_SeqAIJ, 2236 MatGetRowIJ_SeqAIJ, 2237 MatRestoreRowIJ_SeqAIJ, 2238 MatGetColumnIJ_SeqAIJ, 2239 MatRestoreColumnIJ_SeqAIJ, 2240 /*55*/ MatFDColoringCreate_SeqAIJ, 2241 0, 2242 0, 2243 MatPermute_SeqAIJ, 2244 0, 2245 /*60*/ 0, 2246 MatDestroy_SeqAIJ, 2247 MatView_SeqAIJ, 2248 MatGetPetscMaps_Petsc, 2249 0, 2250 /*65*/ 0, 2251 0, 2252 0, 2253 0, 2254 0, 2255 /*70*/ MatGetRowMax_SeqAIJ, 2256 0, 2257 MatSetColoring_SeqAIJ, 2258 #if defined(PETSC_HAVE_ADIC) 2259 MatSetValuesAdic_SeqAIJ, 2260 #else 2261 0, 2262 #endif 2263 MatSetValuesAdifor_SeqAIJ, 2264 /*75*/ MatFDColoringApply_SeqAIJ, 2265 0, 2266 0, 2267 0, 2268 0, 2269 /*80*/ 0, 2270 0, 2271 0, 2272 0, 2273 MatLoad_SeqAIJ, 2274 /*85*/ MatIsSymmetric_SeqAIJ, 2275 0, 2276 0, 2277 0, 2278 0, 2279 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2280 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2281 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2282 MatPtAP_Basic, 2283 MatPtAPSymbolic_SeqAIJ, 2284 /*95*/ MatPtAPNumeric_SeqAIJ, 2285 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2286 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2287 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2288 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2289 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2290 0, 2291 0, 2292 MatConjugate_SeqAIJ 2293 }; 2294 2295 EXTERN_C_BEGIN 2296 #undef __FUNCT__ 2297 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2298 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2299 { 2300 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2301 PetscInt i,nz,n; 2302 2303 PetscFunctionBegin; 2304 2305 nz = aij->maxnz; 2306 n = mat->n; 2307 for (i=0; i<nz; i++) { 2308 aij->j[i] = indices[i]; 2309 } 2310 aij->nz = nz; 2311 for (i=0; i<n; i++) { 2312 aij->ilen[i] = aij->imax[i]; 2313 } 2314 2315 PetscFunctionReturn(0); 2316 } 2317 EXTERN_C_END 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2321 /*@ 2322 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2323 in the matrix. 2324 2325 Input Parameters: 2326 + mat - the SeqAIJ matrix 2327 - indices - the column indices 2328 2329 Level: advanced 2330 2331 Notes: 2332 This can be called if you have precomputed the nonzero structure of the 2333 matrix and want to provide it to the matrix object to improve the performance 2334 of the MatSetValues() operation. 2335 2336 You MUST have set the correct numbers of nonzeros per row in the call to 2337 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2338 2339 MUST be called before any calls to MatSetValues(); 2340 2341 The indices should start with zero, not one. 2342 2343 @*/ 2344 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2345 { 2346 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2347 2348 PetscFunctionBegin; 2349 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2350 PetscValidPointer(indices,2); 2351 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2352 if (f) { 2353 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2354 } else { 2355 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2356 } 2357 PetscFunctionReturn(0); 2358 } 2359 2360 /* ----------------------------------------------------------------------------------------*/ 2361 2362 EXTERN_C_BEGIN 2363 #undef __FUNCT__ 2364 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2365 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2366 { 2367 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2368 PetscErrorCode ierr; 2369 size_t nz = aij->i[mat->m]; 2370 2371 PetscFunctionBegin; 2372 if (aij->nonew != 1) { 2373 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2374 } 2375 2376 /* allocate space for values if not already there */ 2377 if (!aij->saved_values) { 2378 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2379 } 2380 2381 /* copy values over */ 2382 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2383 PetscFunctionReturn(0); 2384 } 2385 EXTERN_C_END 2386 2387 #undef __FUNCT__ 2388 #define __FUNCT__ "MatStoreValues" 2389 /*@ 2390 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2391 example, reuse of the linear part of a Jacobian, while recomputing the 2392 nonlinear portion. 2393 2394 Collect on Mat 2395 2396 Input Parameters: 2397 . mat - the matrix (currently only AIJ matrices support this option) 2398 2399 Level: advanced 2400 2401 Common Usage, with SNESSolve(): 2402 $ Create Jacobian matrix 2403 $ Set linear terms into matrix 2404 $ Apply boundary conditions to matrix, at this time matrix must have 2405 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2406 $ boundary conditions again will not change the nonzero structure 2407 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2408 $ ierr = MatStoreValues(mat); 2409 $ Call SNESSetJacobian() with matrix 2410 $ In your Jacobian routine 2411 $ ierr = MatRetrieveValues(mat); 2412 $ Set nonlinear terms in matrix 2413 2414 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2415 $ // build linear portion of Jacobian 2416 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2417 $ ierr = MatStoreValues(mat); 2418 $ loop over nonlinear iterations 2419 $ ierr = MatRetrieveValues(mat); 2420 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2421 $ // call MatAssemblyBegin/End() on matrix 2422 $ Solve linear system with Jacobian 2423 $ endloop 2424 2425 Notes: 2426 Matrix must already be assemblied before calling this routine 2427 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2428 calling this routine. 2429 2430 When this is called multiple times it overwrites the previous set of stored values 2431 and does not allocated additional space. 2432 2433 .seealso: MatRetrieveValues() 2434 2435 @*/ 2436 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2437 { 2438 PetscErrorCode ierr,(*f)(Mat); 2439 2440 PetscFunctionBegin; 2441 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2442 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2443 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2444 2445 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2446 if (f) { 2447 ierr = (*f)(mat);CHKERRQ(ierr); 2448 } else { 2449 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2450 } 2451 PetscFunctionReturn(0); 2452 } 2453 2454 EXTERN_C_BEGIN 2455 #undef __FUNCT__ 2456 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2457 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2458 { 2459 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2460 PetscErrorCode ierr; 2461 PetscInt nz = aij->i[mat->m]; 2462 2463 PetscFunctionBegin; 2464 if (aij->nonew != 1) { 2465 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2466 } 2467 if (!aij->saved_values) { 2468 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2469 } 2470 /* copy values over */ 2471 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2472 PetscFunctionReturn(0); 2473 } 2474 EXTERN_C_END 2475 2476 #undef __FUNCT__ 2477 #define __FUNCT__ "MatRetrieveValues" 2478 /*@ 2479 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2480 example, reuse of the linear part of a Jacobian, while recomputing the 2481 nonlinear portion. 2482 2483 Collect on Mat 2484 2485 Input Parameters: 2486 . mat - the matrix (currently on AIJ matrices support this option) 2487 2488 Level: advanced 2489 2490 .seealso: MatStoreValues() 2491 2492 @*/ 2493 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2494 { 2495 PetscErrorCode ierr,(*f)(Mat); 2496 2497 PetscFunctionBegin; 2498 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2499 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2500 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2501 2502 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2503 if (f) { 2504 ierr = (*f)(mat);CHKERRQ(ierr); 2505 } else { 2506 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2507 } 2508 PetscFunctionReturn(0); 2509 } 2510 2511 2512 /* --------------------------------------------------------------------------------*/ 2513 #undef __FUNCT__ 2514 #define __FUNCT__ "MatCreateSeqAIJ" 2515 /*@C 2516 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2517 (the default parallel PETSc format). For good matrix assembly performance 2518 the user should preallocate the matrix storage by setting the parameter nz 2519 (or the array nnz). By setting these parameters accurately, performance 2520 during matrix assembly can be increased by more than a factor of 50. 2521 2522 Collective on MPI_Comm 2523 2524 Input Parameters: 2525 + comm - MPI communicator, set to PETSC_COMM_SELF 2526 . m - number of rows 2527 . n - number of columns 2528 . nz - number of nonzeros per row (same for all rows) 2529 - nnz - array containing the number of nonzeros in the various rows 2530 (possibly different for each row) or PETSC_NULL 2531 2532 Output Parameter: 2533 . A - the matrix 2534 2535 Notes: 2536 If nnz is given then nz is ignored 2537 2538 The AIJ format (also called the Yale sparse matrix format or 2539 compressed row storage), is fully compatible with standard Fortran 77 2540 storage. That is, the stored row and column indices can begin at 2541 either one (as in Fortran) or zero. See the users' manual for details. 2542 2543 Specify the preallocated storage with either nz or nnz (not both). 2544 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2545 allocation. For large problems you MUST preallocate memory or you 2546 will get TERRIBLE performance, see the users' manual chapter on matrices. 2547 2548 By default, this format uses inodes (identical nodes) when possible, to 2549 improve numerical efficiency of matrix-vector products and solves. We 2550 search for consecutive rows with the same nonzero structure, thereby 2551 reusing matrix information to achieve increased efficiency. 2552 2553 Options Database Keys: 2554 + -mat_no_inode - Do not use inodes 2555 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2556 - -mat_aij_oneindex - Internally use indexing starting at 1 2557 rather than 0. Note that when calling MatSetValues(), 2558 the user still MUST index entries starting at 0! 2559 2560 Level: intermediate 2561 2562 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2563 2564 @*/ 2565 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2566 { 2567 PetscErrorCode ierr; 2568 2569 PetscFunctionBegin; 2570 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2571 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2572 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2573 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2574 PetscFunctionReturn(0); 2575 } 2576 2577 #undef __FUNCT__ 2578 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2579 /*@C 2580 MatSeqAIJSetPreallocation - For good matrix assembly performance 2581 the user should preallocate the matrix storage by setting the parameter nz 2582 (or the array nnz). By setting these parameters accurately, performance 2583 during matrix assembly can be increased by more than a factor of 50. 2584 2585 Collective on MPI_Comm 2586 2587 Input Parameters: 2588 + B - The matrix 2589 . nz - number of nonzeros per row (same for all rows) 2590 - nnz - array containing the number of nonzeros in the various rows 2591 (possibly different for each row) or PETSC_NULL 2592 2593 Notes: 2594 If nnz is given then nz is ignored 2595 2596 The AIJ format (also called the Yale sparse matrix format or 2597 compressed row storage), is fully compatible with standard Fortran 77 2598 storage. That is, the stored row and column indices can begin at 2599 either one (as in Fortran) or zero. See the users' manual for details. 2600 2601 Specify the preallocated storage with either nz or nnz (not both). 2602 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2603 allocation. For large problems you MUST preallocate memory or you 2604 will get TERRIBLE performance, see the users' manual chapter on matrices. 2605 2606 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2607 entries or columns indices 2608 2609 By default, this format uses inodes (identical nodes) when possible, to 2610 improve numerical efficiency of matrix-vector products and solves. We 2611 search for consecutive rows with the same nonzero structure, thereby 2612 reusing matrix information to achieve increased efficiency. 2613 2614 Options Database Keys: 2615 + -mat_no_inode - Do not use inodes 2616 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2617 - -mat_aij_oneindex - Internally use indexing starting at 1 2618 rather than 0. Note that when calling MatSetValues(), 2619 the user still MUST index entries starting at 0! 2620 2621 Level: intermediate 2622 2623 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2624 2625 @*/ 2626 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2627 { 2628 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2629 2630 PetscFunctionBegin; 2631 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2632 if (f) { 2633 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2634 } 2635 PetscFunctionReturn(0); 2636 } 2637 2638 EXTERN_C_BEGIN 2639 #undef __FUNCT__ 2640 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2641 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2642 { 2643 Mat_SeqAIJ *b; 2644 PetscTruth skipallocation = PETSC_FALSE; 2645 PetscErrorCode ierr; 2646 PetscInt i; 2647 2648 PetscFunctionBegin; 2649 2650 if (nz == MAT_SKIP_ALLOCATION) { 2651 skipallocation = PETSC_TRUE; 2652 nz = 0; 2653 } 2654 2655 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2656 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2657 if (nnz) { 2658 for (i=0; i<B->m; i++) { 2659 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2660 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2661 } 2662 } 2663 2664 B->preallocated = PETSC_TRUE; 2665 b = (Mat_SeqAIJ*)B->data; 2666 2667 if (!skipallocation) { 2668 ierr = PetscMalloc2(B->m,PetscInt,&b->imax,B->m,PetscInt,&b->ilen);CHKERRQ(ierr); 2669 if (!nnz) { 2670 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2671 else if (nz <= 0) nz = 1; 2672 for (i=0; i<B->m; i++) b->imax[i] = nz; 2673 nz = nz*B->m; 2674 } else { 2675 nz = 0; 2676 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2677 } 2678 2679 /* b->ilen will count nonzeros in each row so far. */ 2680 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2681 2682 /* allocate the matrix space */ 2683 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->m+1,PetscInt,&b->i);CHKERRQ(ierr); 2684 b->i[0] = 0; 2685 for (i=1; i<B->m+1; i++) { 2686 b->i[i] = b->i[i-1] + b->imax[i-1]; 2687 } 2688 b->singlemalloc = PETSC_TRUE; 2689 b->freedata = PETSC_TRUE; 2690 } else { 2691 b->freedata = PETSC_FALSE; 2692 } 2693 2694 b->nz = 0; 2695 b->maxnz = nz; 2696 B->info.nz_unneeded = (double)b->maxnz; 2697 PetscFunctionReturn(0); 2698 } 2699 EXTERN_C_END 2700 2701 #undef __FUNCT__ 2702 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2703 /*@C 2704 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2705 2706 Input Parameters: 2707 + B - the matrix 2708 . i - the indices into j for the start of each row (starts with zero) 2709 . j - the column indices for each row (starts with zero) these must be sorted for each row 2710 - v - optional values in the matrix 2711 2712 Contributed by: Lisandro Dalchin 2713 2714 Level: developer 2715 2716 .keywords: matrix, aij, compressed row, sparse, sequential 2717 2718 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2719 @*/ 2720 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2721 { 2722 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2723 PetscErrorCode ierr; 2724 2725 PetscFunctionBegin; 2726 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2727 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2728 if (f) { 2729 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2730 } 2731 PetscFunctionReturn(0); 2732 } 2733 2734 EXTERN_C_BEGIN 2735 #undef __FUNCT__ 2736 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2737 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt I[],const PetscInt J[],const PetscScalar v[]) 2738 { 2739 PetscInt i; 2740 PetscInt m,n; 2741 PetscInt nz; 2742 PetscInt *nnz, nz_max = 0; 2743 PetscScalar *values; 2744 PetscErrorCode ierr; 2745 2746 PetscFunctionBegin; 2747 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2748 2749 if (I[0]) { 2750 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "I[0] must be 0 it is %D", I[0]); 2751 } 2752 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2753 for(i = 0; i < m; i++) { 2754 nz = I[i+1]- I[i]; 2755 nz_max = PetscMax(nz_max, nz); 2756 if (nz < 0) { 2757 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2758 } 2759 nnz[i] = nz; 2760 } 2761 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2762 ierr = PetscFree(nnz);CHKERRQ(ierr); 2763 2764 if (v) { 2765 values = (PetscScalar*) v; 2766 } else { 2767 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2768 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2769 } 2770 2771 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 2772 2773 for(i = 0; i < m; i++) { 2774 nz = I[i+1] - I[i]; 2775 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+I[i], values + (v ? I[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2776 } 2777 2778 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2779 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2780 ierr = MatSetOption(B,MAT_COLUMNS_UNSORTED);CHKERRQ(ierr); 2781 2782 if (!v) { 2783 ierr = PetscFree(values);CHKERRQ(ierr); 2784 } 2785 PetscFunctionReturn(0); 2786 } 2787 EXTERN_C_END 2788 2789 /*MC 2790 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2791 based on compressed sparse row format. 2792 2793 Options Database Keys: 2794 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2795 2796 Level: beginner 2797 2798 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2799 M*/ 2800 2801 EXTERN_C_BEGIN 2802 #undef __FUNCT__ 2803 #define __FUNCT__ "MatCreate_SeqAIJ" 2804 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 2805 { 2806 Mat_SeqAIJ *b; 2807 PetscErrorCode ierr; 2808 PetscMPIInt size; 2809 2810 PetscFunctionBegin; 2811 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2812 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2813 2814 B->m = B->M = PetscMax(B->m,B->M); 2815 B->n = B->N = PetscMax(B->n,B->N); 2816 2817 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2818 B->data = (void*)b; 2819 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2820 B->factor = 0; 2821 B->mapping = 0; 2822 b->row = 0; 2823 b->col = 0; 2824 b->icol = 0; 2825 b->reallocs = 0; 2826 2827 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2828 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2829 2830 b->sorted = PETSC_FALSE; 2831 b->ignorezeroentries = PETSC_FALSE; 2832 b->roworiented = PETSC_TRUE; 2833 b->nonew = 0; 2834 b->diag = 0; 2835 b->solve_work = 0; 2836 B->spptr = 0; 2837 b->saved_values = 0; 2838 b->idiag = 0; 2839 b->ssor = 0; 2840 b->keepzeroedrows = PETSC_FALSE; 2841 b->xtoy = 0; 2842 b->XtoY = 0; 2843 b->compressedrow.use = PETSC_FALSE; 2844 b->compressedrow.nrows = B->m; 2845 b->compressedrow.i = PETSC_NULL; 2846 b->compressedrow.rindex = PETSC_NULL; 2847 b->compressedrow.checked = PETSC_FALSE; 2848 B->same_nonzero = PETSC_FALSE; 2849 2850 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2851 2852 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2853 "MatSeqAIJSetColumnIndices_SeqAIJ", 2854 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2855 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2856 "MatStoreValues_SeqAIJ", 2857 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2858 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2859 "MatRetrieveValues_SeqAIJ", 2860 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2861 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2862 "MatConvert_SeqAIJ_SeqSBAIJ", 2863 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2864 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2865 "MatConvert_SeqAIJ_SeqBAIJ", 2866 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2867 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2868 "MatIsTranspose_SeqAIJ", 2869 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2870 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2871 "MatSeqAIJSetPreallocation_SeqAIJ", 2872 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2873 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 2874 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 2875 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 2876 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2877 "MatReorderForNonzeroDiagonal_SeqAIJ", 2878 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2879 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 2880 PetscFunctionReturn(0); 2881 } 2882 EXTERN_C_END 2883 2884 #undef __FUNCT__ 2885 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2886 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2887 { 2888 Mat C; 2889 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2890 PetscErrorCode ierr; 2891 PetscInt i,m = A->m; 2892 2893 PetscFunctionBegin; 2894 *B = 0; 2895 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 2896 ierr = MatSetSizes(C,A->m,A->n,A->m,A->n);CHKERRQ(ierr); 2897 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2898 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2899 2900 c = (Mat_SeqAIJ*)C->data; 2901 2902 C->factor = A->factor; 2903 2904 c->row = 0; 2905 c->col = 0; 2906 c->icol = 0; 2907 c->reallocs = 0; 2908 2909 C->assembled = PETSC_TRUE; 2910 2911 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 2912 for (i=0; i<m; i++) { 2913 c->imax[i] = a->imax[i]; 2914 c->ilen[i] = a->ilen[i]; 2915 } 2916 2917 /* allocate the matrix space */ 2918 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 2919 c->singlemalloc = PETSC_TRUE; 2920 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2921 if (m > 0) { 2922 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2923 if (cpvalues == MAT_COPY_VALUES) { 2924 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2925 } else { 2926 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2927 } 2928 } 2929 2930 c->sorted = a->sorted; 2931 c->ignorezeroentries = a->ignorezeroentries; 2932 c->roworiented = a->roworiented; 2933 c->nonew = a->nonew; 2934 if (a->diag) { 2935 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2936 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2937 for (i=0; i<m; i++) { 2938 c->diag[i] = a->diag[i]; 2939 } 2940 } else c->diag = 0; 2941 c->solve_work = 0; 2942 c->saved_values = 0; 2943 c->idiag = 0; 2944 c->ssor = 0; 2945 c->keepzeroedrows = a->keepzeroedrows; 2946 c->freedata = PETSC_TRUE; 2947 c->xtoy = 0; 2948 c->XtoY = 0; 2949 2950 c->nz = a->nz; 2951 c->maxnz = a->maxnz; 2952 C->preallocated = PETSC_TRUE; 2953 2954 c->compressedrow.use = a->compressedrow.use; 2955 c->compressedrow.nrows = a->compressedrow.nrows; 2956 c->compressedrow.checked = a->compressedrow.checked; 2957 if ( a->compressedrow.checked && a->compressedrow.use){ 2958 i = a->compressedrow.nrows; 2959 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 2960 c->compressedrow.rindex = c->compressedrow.i + i + 1; 2961 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 2962 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 2963 } else { 2964 c->compressedrow.use = PETSC_FALSE; 2965 c->compressedrow.i = PETSC_NULL; 2966 c->compressedrow.rindex = PETSC_NULL; 2967 } 2968 C->same_nonzero = A->same_nonzero; 2969 2970 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 2971 2972 *B = C; 2973 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2974 PetscFunctionReturn(0); 2975 } 2976 2977 #undef __FUNCT__ 2978 #define __FUNCT__ "MatLoad_SeqAIJ" 2979 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 2980 { 2981 Mat_SeqAIJ *a; 2982 Mat B; 2983 PetscErrorCode ierr; 2984 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 2985 int fd; 2986 PetscMPIInt size; 2987 MPI_Comm comm; 2988 2989 PetscFunctionBegin; 2990 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2991 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2992 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2993 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2994 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2995 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2996 M = header[1]; N = header[2]; nz = header[3]; 2997 2998 if (nz < 0) { 2999 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3000 } 3001 3002 /* read in row lengths */ 3003 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3004 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3005 3006 /* check if sum of rowlengths is same as nz */ 3007 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3008 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3009 3010 /* create our matrix */ 3011 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3012 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3013 ierr = MatSetType(B,type);CHKERRQ(ierr); 3014 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3015 a = (Mat_SeqAIJ*)B->data; 3016 3017 /* read in column indices and adjust for Fortran indexing*/ 3018 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3019 3020 /* read in nonzero values */ 3021 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3022 3023 /* set matrix "i" values */ 3024 a->i[0] = 0; 3025 for (i=1; i<= M; i++) { 3026 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3027 a->ilen[i-1] = rowlengths[i-1]; 3028 } 3029 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3030 3031 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3032 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3033 *A = B; 3034 PetscFunctionReturn(0); 3035 } 3036 3037 #undef __FUNCT__ 3038 #define __FUNCT__ "MatEqual_SeqAIJ" 3039 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3040 { 3041 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3042 PetscErrorCode ierr; 3043 3044 PetscFunctionBegin; 3045 /* If the matrix dimensions are not equal,or no of nonzeros */ 3046 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 3047 *flg = PETSC_FALSE; 3048 PetscFunctionReturn(0); 3049 } 3050 3051 /* if the a->i are the same */ 3052 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3053 if (!*flg) PetscFunctionReturn(0); 3054 3055 /* if a->j are the same */ 3056 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3057 if (!*flg) PetscFunctionReturn(0); 3058 3059 /* if a->a are the same */ 3060 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3061 3062 PetscFunctionReturn(0); 3063 3064 } 3065 3066 #undef __FUNCT__ 3067 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3068 /*@ 3069 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3070 provided by the user. 3071 3072 Coolective on MPI_Comm 3073 3074 Input Parameters: 3075 + comm - must be an MPI communicator of size 1 3076 . m - number of rows 3077 . n - number of columns 3078 . i - row indices 3079 . j - column indices 3080 - a - matrix values 3081 3082 Output Parameter: 3083 . mat - the matrix 3084 3085 Level: intermediate 3086 3087 Notes: 3088 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3089 once the matrix is destroyed 3090 3091 You cannot set new nonzero locations into this matrix, that will generate an error. 3092 3093 The i and j indices are 0 based 3094 3095 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 3096 3097 @*/ 3098 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3099 { 3100 PetscErrorCode ierr; 3101 PetscInt ii; 3102 Mat_SeqAIJ *aij; 3103 3104 PetscFunctionBegin; 3105 if (i[0]) { 3106 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3107 } 3108 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3109 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3110 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3111 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3112 aij = (Mat_SeqAIJ*)(*mat)->data; 3113 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3114 3115 aij->i = i; 3116 aij->j = j; 3117 aij->a = a; 3118 aij->singlemalloc = PETSC_FALSE; 3119 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3120 aij->freedata = PETSC_FALSE; 3121 3122 for (ii=0; ii<m; ii++) { 3123 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3124 #if defined(PETSC_USE_DEBUG) 3125 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3126 #endif 3127 } 3128 #if defined(PETSC_USE_DEBUG) 3129 for (ii=0; ii<aij->i[m]; ii++) { 3130 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3131 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3132 } 3133 #endif 3134 3135 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3136 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3137 PetscFunctionReturn(0); 3138 } 3139 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3142 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3143 { 3144 PetscErrorCode ierr; 3145 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3146 3147 PetscFunctionBegin; 3148 if (coloring->ctype == IS_COLORING_LOCAL) { 3149 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3150 a->coloring = coloring; 3151 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3152 PetscInt i,*larray; 3153 ISColoring ocoloring; 3154 ISColoringValue *colors; 3155 3156 /* set coloring for diagonal portion */ 3157 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3158 for (i=0; i<A->n; i++) { 3159 larray[i] = i; 3160 } 3161 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3162 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3163 for (i=0; i<A->n; i++) { 3164 colors[i] = coloring->colors[larray[i]]; 3165 } 3166 ierr = PetscFree(larray);CHKERRQ(ierr); 3167 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3168 a->coloring = ocoloring; 3169 } 3170 PetscFunctionReturn(0); 3171 } 3172 3173 #if defined(PETSC_HAVE_ADIC) 3174 EXTERN_C_BEGIN 3175 #include "adic/ad_utils.h" 3176 EXTERN_C_END 3177 3178 #undef __FUNCT__ 3179 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3180 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3181 { 3182 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3183 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3184 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3185 ISColoringValue *color; 3186 3187 PetscFunctionBegin; 3188 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3189 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3190 color = a->coloring->colors; 3191 /* loop over rows */ 3192 for (i=0; i<m; i++) { 3193 nz = ii[i+1] - ii[i]; 3194 /* loop over columns putting computed value into matrix */ 3195 for (j=0; j<nz; j++) { 3196 *v++ = values[color[*jj++]]; 3197 } 3198 values += nlen; /* jump to next row of derivatives */ 3199 } 3200 PetscFunctionReturn(0); 3201 } 3202 #endif 3203 3204 #undef __FUNCT__ 3205 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3206 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3207 { 3208 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3209 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3210 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3211 ISColoringValue *color; 3212 3213 PetscFunctionBegin; 3214 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3215 color = a->coloring->colors; 3216 /* loop over rows */ 3217 for (i=0; i<m; i++) { 3218 nz = ii[i+1] - ii[i]; 3219 /* loop over columns putting computed value into matrix */ 3220 for (j=0; j<nz; j++) { 3221 *v++ = values[color[*jj++]]; 3222 } 3223 values += nl; /* jump to next row of derivatives */ 3224 } 3225 PetscFunctionReturn(0); 3226 } 3227 3228