Actual source code: sbaijfact.c
petsc-3.3-p7 2013-05-11
2: #include <../src/mat/impls/baij/seq/baij.h>
3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
4: #include <../src/mat/blockinvert.h>
5: #include <petscis.h>
7: /*
8: input:
9: F -- numeric factor
10: output:
11: nneg, nzero, npos: matrix inertia
12: */
16: PetscErrorCode MatGetInertia_SeqSBAIJ(Mat F,PetscInt *nneig,PetscInt *nzero,PetscInt *npos)
17: {
18: Mat_SeqSBAIJ *fact_ptr = (Mat_SeqSBAIJ*)F->data;
19: MatScalar *dd = fact_ptr->a;
20: PetscInt mbs=fact_ptr->mbs,bs=F->rmap->bs,i,nneig_tmp,npos_tmp,*fi = fact_ptr->diag;
23: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for bs: %D >1 yet",bs);
24: nneig_tmp = 0; npos_tmp = 0;
25: for (i=0; i<mbs; i++){
26: if (PetscRealPart(dd[*fi]) > 0.0){
27: npos_tmp++;
28: } else if (PetscRealPart(dd[*fi]) < 0.0){
29: nneig_tmp++;
30: }
31: fi++;
32: }
33: if (nneig) *nneig = nneig_tmp;
34: if (npos) *npos = npos_tmp;
35: if (nzero) *nzero = mbs - nneig_tmp - npos_tmp;
37: return(0);
38: }
40: /*
41: Symbolic U^T*D*U factorization for SBAIJ format. Modified from SSF of YSMP.
42: Use Modified Sparse Row (MSR) storage for u and ju. See page 85, "Iterative Methods ..." by Saad.
43: */
46: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(Mat F,Mat A,IS perm,const MatFactorInfo *info)
47: {
48: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b;
50: const PetscInt *rip,*ai,*aj;
51: PetscInt i,mbs = a->mbs,*jutmp,bs = A->rmap->bs,bs2=a->bs2;
52: PetscInt m,reallocs = 0,prow;
53: PetscInt *jl,*q,jmin,jmax,juidx,nzk,qm,*iu,*ju,k,j,vj,umax,maxadd;
54: PetscReal f = info->fill;
55: PetscBool perm_identity;
58: /* check whether perm is the identity mapping */
59: ISIdentity(perm,&perm_identity);
60: ISGetIndices(perm,&rip);
61:
62: if (perm_identity){ /* without permutation */
63: a->permute = PETSC_FALSE;
64: ai = a->i; aj = a->j;
65: } else { /* non-trivial permutation */
66: a->permute = PETSC_TRUE;
67: MatReorderingSeqSBAIJ(A,perm);
68: ai = a->inew; aj = a->jnew;
69: }
70:
71: /* initialization */
72: PetscMalloc((mbs+1)*sizeof(PetscInt),&iu);
73: umax = (PetscInt)(f*ai[mbs] + 1); umax += mbs + 1;
74: PetscMalloc(umax*sizeof(PetscInt),&ju);
75: iu[0] = mbs+1;
76: juidx = mbs + 1; /* index for ju */
77: /* jl linked list for pivot row -- linked list for col index */
78: PetscMalloc2(mbs,PetscInt,&jl,mbs,PetscInt,&q);
79: for (i=0; i<mbs; i++){
80: jl[i] = mbs;
81: q[i] = 0;
82: }
84: /* for each row k */
85: for (k=0; k<mbs; k++){
86: for (i=0; i<mbs; i++) q[i] = 0; /* to be removed! */
87: nzk = 0; /* num. of nz blocks in k-th block row with diagonal block excluded */
88: q[k] = mbs;
89: /* initialize nonzero structure of k-th row to row rip[k] of A */
90: jmin = ai[rip[k]] +1; /* exclude diag[k] */
91: jmax = ai[rip[k]+1];
92: for (j=jmin; j<jmax; j++){
93: vj = rip[aj[j]]; /* col. value */
94: if(vj > k){
95: qm = k;
96: do {
97: m = qm; qm = q[m];
98: } while(qm < vj);
99: if (qm == vj) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Duplicate entry in A\n");
100: nzk++;
101: q[m] = vj;
102: q[vj] = qm;
103: } /* if(vj > k) */
104: } /* for (j=jmin; j<jmax; j++) */
106: /* modify nonzero structure of k-th row by computing fill-in
107: for each row i to be merged in */
108: prow = k;
109: prow = jl[prow]; /* next pivot row (== mbs for symbolic factorization) */
110:
111: while (prow < k){
112: /* merge row prow into k-th row */
113: jmin = iu[prow] + 1; jmax = iu[prow+1];
114: qm = k;
115: for (j=jmin; j<jmax; j++){
116: vj = ju[j];
117: do {
118: m = qm; qm = q[m];
119: } while (qm < vj);
120: if (qm != vj){
121: nzk++; q[m] = vj; q[vj] = qm; qm = vj;
122: }
123: }
124: prow = jl[prow]; /* next pivot row */
125: }
126:
127: /* add k to row list for first nonzero element in k-th row */
128: if (nzk > 0){
129: i = q[k]; /* col value of first nonzero element in U(k, k+1:mbs-1) */
130: jl[k] = jl[i]; jl[i] = k;
131: }
132: iu[k+1] = iu[k] + nzk;
134: /* allocate more space to ju if needed */
135: if (iu[k+1] > umax) {
136: /* estimate how much additional space we will need */
137: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
138: /* just double the memory each time */
139: maxadd = umax;
140: if (maxadd < nzk) maxadd = (mbs-k)*(nzk+1)/2;
141: umax += maxadd;
143: /* allocate a longer ju */
144: PetscMalloc(umax*sizeof(PetscInt),&jutmp);
145: PetscMemcpy(jutmp,ju,iu[k]*sizeof(PetscInt));
146: PetscFree(ju);
147: ju = jutmp;
148: reallocs++; /* count how many times we realloc */
149: }
151: /* save nonzero structure of k-th row in ju */
152: i=k;
153: while (nzk --) {
154: i = q[i];
155: ju[juidx++] = i;
156: }
157: }
159: #if defined(PETSC_USE_INFO)
160: if (ai[mbs] != 0) {
161: PetscReal af = ((PetscReal)iu[mbs])/((PetscReal)ai[mbs]);
162: PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);
163: PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);
164: PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);
165: PetscInfo(A,"for best performance.\n");
166: } else {
167: PetscInfo(A,"Empty matrix.\n");
168: }
169: #endif
171: ISRestoreIndices(perm,&rip);
172: PetscFree2(jl,q);
174: /* put together the new matrix */
175: MatSeqSBAIJSetPreallocation_SeqSBAIJ(F,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);
177: /* PetscLogObjectParent(B,iperm); */
178: b = (Mat_SeqSBAIJ*)(F)->data;
179: b->singlemalloc = PETSC_FALSE;
180: b->free_a = PETSC_TRUE;
181: b->free_ij = PETSC_TRUE;
182: PetscMalloc((iu[mbs]+1)*sizeof(MatScalar)*bs2,&b->a);
183: b->j = ju;
184: b->i = iu;
185: b->diag = 0;
186: b->ilen = 0;
187: b->imax = 0;
188: b->row = perm;
189: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
190: PetscObjectReference((PetscObject)perm);
191: b->icol = perm;
192: PetscObjectReference((PetscObject)perm);
193: PetscMalloc((bs*mbs+bs)*sizeof(PetscScalar),&b->solve_work);
194: /* In b structure: Free imax, ilen, old a, old j.
195: Allocate idnew, solve_work, new a, new j */
196: PetscLogObjectMemory(F,(iu[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
197: b->maxnz = b->nz = iu[mbs];
198:
199: (F)->info.factor_mallocs = reallocs;
200: (F)->info.fill_ratio_given = f;
201: if (ai[mbs] != 0) {
202: (F)->info.fill_ratio_needed = ((PetscReal)iu[mbs])/((PetscReal)ai[mbs]);
203: } else {
204: (F)->info.fill_ratio_needed = 0.0;
205: }
206: MatSeqSBAIJSetNumericFactorization_inplace(F,perm_identity);
207: return(0);
208: }
209: /*
210: Symbolic U^T*D*U factorization for SBAIJ format.
211: See MatICCFactorSymbolic_SeqAIJ() for description of its data structure.
212: */
213: #include <petscbt.h>
214: #include <../src/mat/utils/freespace.h>
217: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
218: {
219: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
220: Mat_SeqSBAIJ *b;
221: PetscErrorCode ierr;
222: PetscBool perm_identity,missing;
223: PetscReal fill = info->fill;
224: const PetscInt *rip,*ai=a->i,*aj=a->j;
225: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,reallocs=0,prow,d;
226: PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
227: PetscInt nlnk,*lnk,ncols,*cols,*uj,**ui_ptr,*uj_ptr,*udiag;
228: PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
229: PetscBT lnkbt;
232: if (bs > 1){
233: MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(fact,A,perm,info);
234: return(0);
235: }
236: if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap->n,A->cmap->n);
237: MatMissingDiagonal(A,&missing,&d);
238: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
240: /* check whether perm is the identity mapping */
241: ISIdentity(perm,&perm_identity);
242: if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
243: a->permute = PETSC_FALSE;
244: ISGetIndices(perm,&rip);
246: /* initialization */
247: PetscMalloc((mbs+1)*sizeof(PetscInt),&ui);
248: PetscMalloc((mbs+1)*sizeof(PetscInt),&udiag);
249: ui[0] = 0;
251: /* jl: linked list for storing indices of the pivot rows
252: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
253: PetscMalloc4(mbs,PetscInt*,&ui_ptr,mbs,PetscInt,&il,mbs,PetscInt,&jl,mbs,PetscInt,&cols);
254: for (i=0; i<mbs; i++){
255: jl[i] = mbs; il[i] = 0;
256: }
258: /* create and initialize a linked list for storing column indices of the active row k */
259: nlnk = mbs + 1;
260: PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);
262: /* initial FreeSpace size is fill*(ai[mbs]+1) */
263: PetscFreeSpaceGet((PetscInt)(fill*(ai[mbs]+1)),&free_space);
264: current_space = free_space;
266: for (k=0; k<mbs; k++){ /* for each active row k */
267: /* initialize lnk by the column indices of row rip[k] of A */
268: nzk = 0;
269: ncols = ai[k+1] - ai[k];
270: if (!ncols) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row %D in matrix ",k);
271: for (j=0; j<ncols; j++){
272: i = *(aj + ai[k] + j);
273: cols[j] = i;
274: }
275: PetscLLAdd(ncols,cols,mbs,nlnk,lnk,lnkbt);
276: nzk += nlnk;
278: /* update lnk by computing fill-in for each pivot row to be merged in */
279: prow = jl[k]; /* 1st pivot row */
280:
281: while (prow < k){
282: nextprow = jl[prow];
283: /* merge prow into k-th row */
284: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
285: jmax = ui[prow+1];
286: ncols = jmax-jmin;
287: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
288: PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
289: nzk += nlnk;
291: /* update il and jl for prow */
292: if (jmin < jmax){
293: il[prow] = jmin;
294: j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
295: }
296: prow = nextprow;
297: }
299: /* if free space is not available, make more free space */
300: if (current_space->local_remaining<nzk) {
301: i = mbs - k + 1; /* num of unfactored rows */
302: i *= PetscMin(nzk, i-1); /* i*nzk, i*(i-1): estimated and max additional space needed */
303: PetscFreeSpaceGet(i,¤t_space);
304: reallocs++;
305: }
307: /* copy data into free space, then initialize lnk */
308: PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);
310: /* add the k-th row into il and jl */
311: if (nzk > 1){
312: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
313: jl[k] = jl[i]; jl[i] = k;
314: il[k] = ui[k] + 1;
315: }
316: ui_ptr[k] = current_space->array;
317: current_space->array += nzk;
318: current_space->local_used += nzk;
319: current_space->local_remaining -= nzk;
321: ui[k+1] = ui[k] + nzk;
322: }
324: ISRestoreIndices(perm,&rip);
325: PetscFree4(ui_ptr,il,jl,cols);
327: /* destroy list of free space and other temporary array(s) */
328: PetscMalloc((ui[mbs]+1)*sizeof(PetscInt),&uj);
329: PetscFreeSpaceContiguous_Cholesky(&free_space,uj,mbs,ui,udiag); /* store matrix factor */
330: PetscLLDestroy(lnk,lnkbt);
332: /* put together the new matrix in MATSEQSBAIJ format */
333: MatSeqSBAIJSetPreallocation_SeqSBAIJ(fact,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);
334:
335: b = (Mat_SeqSBAIJ*)fact->data;
336: b->singlemalloc = PETSC_FALSE;
337: b->free_a = PETSC_TRUE;
338: b->free_ij = PETSC_TRUE;
339: PetscMalloc((ui[mbs]+1)*sizeof(MatScalar),&b->a);
340: b->j = uj;
341: b->i = ui;
342: b->diag = udiag;
343: b->free_diag = PETSC_TRUE;
344: b->ilen = 0;
345: b->imax = 0;
346: b->row = perm;
347: b->icol = perm;
348: PetscObjectReference((PetscObject)perm);
349: PetscObjectReference((PetscObject)perm);
350: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
351: PetscMalloc((mbs+1)*sizeof(PetscScalar),&b->solve_work);
352: PetscLogObjectMemory(fact,ui[mbs]*(sizeof(PetscInt)+sizeof(MatScalar)));
353: b->maxnz = b->nz = ui[mbs];
354:
355: fact->info.factor_mallocs = reallocs;
356: fact->info.fill_ratio_given = fill;
357: if (ai[mbs] != 0) {
358: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs])/ai[mbs];
359: } else {
360: fact->info.fill_ratio_needed = 0.0;
361: }
362: #if defined(PETSC_USE_INFO)
363: if (ai[mbs] != 0) {
364: PetscReal af = fact->info.fill_ratio_needed;
365: PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);
366: PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);
367: PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);
368: } else {
369: PetscInfo(A,"Empty matrix.\n");
370: }
371: #endif
372: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering;
373: return(0);
374: }
378: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
379: {
380: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
381: Mat_SeqSBAIJ *b;
382: PetscErrorCode ierr;
383: PetscBool perm_identity,missing;
384: PetscReal fill = info->fill;
385: const PetscInt *rip,*ai,*aj;
386: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,reallocs=0,prow,d;
387: PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
388: PetscInt nlnk,*lnk,ncols,*cols,*uj,**ui_ptr,*uj_ptr;
389: PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
390: PetscBT lnkbt;
393: MatMissingDiagonal(A,&missing,&d);
394: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
396: /*
397: This code originally uses Modified Sparse Row (MSR) storage
398: (see page 85, "Iterative Methods ..." by Saad) for the output matrix B - bad choise!
399: Then it is rewritten so the factor B takes seqsbaij format. However the associated
400: MatCholeskyFactorNumeric_() have not been modified for the cases of bs>1 or !perm_identity,
401: thus the original code in MSR format is still used for these cases.
402: The code below should replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR() whenever
403: MatCholeskyFactorNumeric_() is modified for using sbaij symbolic factor.
404: */
405: if (bs > 1){
406: MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(fact,A,perm,info);
407: return(0);
408: }
410: /* check whether perm is the identity mapping */
411: ISIdentity(perm,&perm_identity);
413: if (perm_identity){
414: a->permute = PETSC_FALSE;
415: ai = a->i; aj = a->j;
416: } else {
417: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
418: /* There are bugs for reordeing. Needs further work.
419: MatReordering for sbaij cannot be efficient. User should use aij formt! */
420: a->permute = PETSC_TRUE;
421: MatReorderingSeqSBAIJ(A,perm);
422: ai = a->inew; aj = a->jnew;
423: }
424: ISGetIndices(perm,&rip);
426: /* initialization */
427: PetscMalloc((mbs+1)*sizeof(PetscInt),&ui);
428: ui[0] = 0;
430: /* jl: linked list for storing indices of the pivot rows
431: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
432: PetscMalloc4(mbs,PetscInt*,&ui_ptr,mbs,PetscInt,&il,mbs,PetscInt,&jl,mbs,PetscInt,&cols);
433: for (i=0; i<mbs; i++){
434: jl[i] = mbs; il[i] = 0;
435: }
437: /* create and initialize a linked list for storing column indices of the active row k */
438: nlnk = mbs + 1;
439: PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);
441: /* initial FreeSpace size is fill*(ai[mbs]+1) */
442: PetscFreeSpaceGet((PetscInt)(fill*(ai[mbs]+1)),&free_space);
443: current_space = free_space;
445: for (k=0; k<mbs; k++){ /* for each active row k */
446: /* initialize lnk by the column indices of row rip[k] of A */
447: nzk = 0;
448: ncols = ai[rip[k]+1] - ai[rip[k]];
449: for (j=0; j<ncols; j++){
450: i = *(aj + ai[rip[k]] + j);
451: cols[j] = rip[i];
452: }
453: PetscLLAdd(ncols,cols,mbs,nlnk,lnk,lnkbt);
454: nzk += nlnk;
456: /* update lnk by computing fill-in for each pivot row to be merged in */
457: prow = jl[k]; /* 1st pivot row */
458:
459: while (prow < k){
460: nextprow = jl[prow];
461: /* merge prow into k-th row */
462: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
463: jmax = ui[prow+1];
464: ncols = jmax-jmin;
465: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
466: PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
467: nzk += nlnk;
469: /* update il and jl for prow */
470: if (jmin < jmax){
471: il[prow] = jmin;
472: j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
473: }
474: prow = nextprow;
475: }
477: /* if free space is not available, make more free space */
478: if (current_space->local_remaining<nzk) {
479: i = mbs - k + 1; /* num of unfactored rows */
480: i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
481: PetscFreeSpaceGet(i,¤t_space);
482: reallocs++;
483: }
485: /* copy data into free space, then initialize lnk */
486: PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);
488: /* add the k-th row into il and jl */
489: if (nzk-1 > 0){
490: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
491: jl[k] = jl[i]; jl[i] = k;
492: il[k] = ui[k] + 1;
493: }
494: ui_ptr[k] = current_space->array;
495: current_space->array += nzk;
496: current_space->local_used += nzk;
497: current_space->local_remaining -= nzk;
499: ui[k+1] = ui[k] + nzk;
500: }
502: ISRestoreIndices(perm,&rip);
503: PetscFree4(ui_ptr,il,jl,cols);
505: /* destroy list of free space and other temporary array(s) */
506: PetscMalloc((ui[mbs]+1)*sizeof(PetscInt),&uj);
507: PetscFreeSpaceContiguous(&free_space,uj);
508: PetscLLDestroy(lnk,lnkbt);
510: /* put together the new matrix in MATSEQSBAIJ format */
511: MatSeqSBAIJSetPreallocation_SeqSBAIJ(fact,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);
512:
513: b = (Mat_SeqSBAIJ*)fact->data;
514: b->singlemalloc = PETSC_FALSE;
515: b->free_a = PETSC_TRUE;
516: b->free_ij = PETSC_TRUE;
517: PetscMalloc((ui[mbs]+1)*sizeof(MatScalar),&b->a);
518: b->j = uj;
519: b->i = ui;
520: b->diag = 0;
521: b->ilen = 0;
522: b->imax = 0;
523: b->row = perm;
524: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
525: PetscObjectReference((PetscObject)perm);
526: b->icol = perm;
527: PetscObjectReference((PetscObject)perm);
528: PetscMalloc((mbs+1)*sizeof(PetscScalar),&b->solve_work);
529: PetscLogObjectMemory(fact,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
530: b->maxnz = b->nz = ui[mbs];
531:
532: fact->info.factor_mallocs = reallocs;
533: fact->info.fill_ratio_given = fill;
534: if (ai[mbs] != 0) {
535: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs])/ai[mbs];
536: } else {
537: fact->info.fill_ratio_needed = 0.0;
538: }
539: #if defined(PETSC_USE_INFO)
540: if (ai[mbs] != 0) {
541: PetscReal af = fact->info.fill_ratio_needed;
542: PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);
543: PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);
544: PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);
545: } else {
546: PetscInfo(A,"Empty matrix.\n");
547: }
548: #endif
549: MatSeqSBAIJSetNumericFactorization_inplace(fact,perm_identity);
550: return(0);
551: }
555: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
556: {
557: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ *)C->data;
558: IS perm = b->row;
560: const PetscInt *ai,*aj,*perm_ptr,mbs=a->mbs,*bi=b->i,*bj=b->j;
561: PetscInt i,j;
562: PetscInt *a2anew,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
563: PetscInt bs=A->rmap->bs,bs2 = a->bs2,bslog = 0;
564: MatScalar *ba = b->a,*aa,*ap,*dk,*uik;
565: MatScalar *u,*diag,*rtmp,*rtmp_ptr;
566: MatScalar *work;
567: PetscInt *pivots;
570: /* initialization */
571: PetscMalloc(bs2*mbs*sizeof(MatScalar),&rtmp);
572: PetscMemzero(rtmp,bs2*mbs*sizeof(MatScalar));
573: PetscMalloc2(mbs,PetscInt,&il,mbs,PetscInt,&jl);
574: for (i=0; i<mbs; i++) {
575: jl[i] = mbs; il[0] = 0;
576: }
577: PetscMalloc3(bs2,MatScalar,&dk,bs2,MatScalar,&uik,bs,MatScalar,&work);
578: PetscMalloc(bs*sizeof(PetscInt),&pivots);
579:
580: ISGetIndices(perm,&perm_ptr);
581:
582: /* check permutation */
583: if (!a->permute){
584: ai = a->i; aj = a->j; aa = a->a;
585: } else {
586: ai = a->inew; aj = a->jnew;
587: PetscMalloc(bs2*ai[mbs]*sizeof(MatScalar),&aa);
588: PetscMemcpy(aa,a->a,bs2*ai[mbs]*sizeof(MatScalar));
589: PetscMalloc(ai[mbs]*sizeof(PetscInt),&a2anew);
590: PetscMemcpy(a2anew,a->a2anew,(ai[mbs])*sizeof(PetscInt));
592: /* flops in while loop */
593: bslog = 2*bs*bs2;
595: for (i=0; i<mbs; i++){
596: jmin = ai[i]; jmax = ai[i+1];
597: for (j=jmin; j<jmax; j++){
598: while (a2anew[j] != j){
599: k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
600: for (k1=0; k1<bs2; k1++){
601: dk[k1] = aa[k*bs2+k1];
602: aa[k*bs2+k1] = aa[j*bs2+k1];
603: aa[j*bs2+k1] = dk[k1];
604: }
605: }
606: /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
607: if (i > aj[j]){
608: /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
609: ap = aa + j*bs2; /* ptr to the beginning of j-th block of aa */
610: for (k=0; k<bs2; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
611: for (k=0; k<bs; k++){ /* j-th block of aa <- dk^T */
612: for (k1=0; k1<bs; k1++) *ap++ = dk[k + bs*k1];
613: }
614: }
615: }
616: }
617: PetscFree(a2anew);
618: }
619:
620: /* for each row k */
621: for (k = 0; k<mbs; k++){
623: /*initialize k-th row with elements nonzero in row perm(k) of A */
624: jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
625:
626: ap = aa + jmin*bs2;
627: for (j = jmin; j < jmax; j++){
628: vj = perm_ptr[aj[j]]; /* block col. index */
629: rtmp_ptr = rtmp + vj*bs2;
630: for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
631: }
633: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
634: PetscMemcpy(dk,rtmp+k*bs2,bs2*sizeof(MatScalar));
635: i = jl[k]; /* first row to be added to k_th row */
637: while (i < k){
638: nexti = jl[i]; /* next row to be added to k_th row */
640: /* compute multiplier */
641: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
643: /* uik = -inv(Di)*U_bar(i,k) */
644: diag = ba + i*bs2;
645: u = ba + ili*bs2;
646: PetscMemzero(uik,bs2*sizeof(MatScalar));
647: PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);
648:
649: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
650: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
651: PetscLogFlops(bslog*2.0);
652:
653: /* update -U(i,k) */
654: PetscMemcpy(ba+ili*bs2,uik,bs2*sizeof(MatScalar));
656: /* add multiple of row i to k-th row ... */
657: jmin = ili + 1; jmax = bi[i+1];
658: if (jmin < jmax){
659: for (j=jmin; j<jmax; j++) {
660: /* rtmp += -U(i,k)^T * U_bar(i,j) */
661: rtmp_ptr = rtmp + bj[j]*bs2;
662: u = ba + j*bs2;
663: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
664: }
665: PetscLogFlops(bslog*(jmax-jmin));
666:
667: /* ... add i to row list for next nonzero entry */
668: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
669: j = bj[jmin];
670: jl[i] = jl[j]; jl[j] = i; /* update jl */
671: }
672: i = nexti;
673: }
675: /* save nonzero entries in k-th row of U ... */
677: /* invert diagonal block */
678: diag = ba+k*bs2;
679: PetscMemcpy(diag,dk,bs2*sizeof(MatScalar));
680: PetscKernel_A_gets_inverse_A(bs,diag,pivots,work);
681:
682: jmin = bi[k]; jmax = bi[k+1];
683: if (jmin < jmax) {
684: for (j=jmin; j<jmax; j++){
685: vj = bj[j]; /* block col. index of U */
686: u = ba + j*bs2;
687: rtmp_ptr = rtmp + vj*bs2;
688: for (k1=0; k1<bs2; k1++){
689: *u++ = *rtmp_ptr;
690: *rtmp_ptr++ = 0.0;
691: }
692: }
693:
694: /* ... add k to row list for first nonzero entry in k-th row */
695: il[k] = jmin;
696: i = bj[jmin];
697: jl[k] = jl[i]; jl[i] = k;
698: }
699: }
701: PetscFree(rtmp);
702: PetscFree2(il,jl);
703: PetscFree3(dk,uik,work);
704: PetscFree(pivots);
705: if (a->permute){
706: PetscFree(aa);
707: }
709: ISRestoreIndices(perm,&perm_ptr);
710: C->ops->solve = MatSolve_SeqSBAIJ_N_inplace;
711: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_inplace;
712: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_inplace;
713: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_inplace;
715: C->assembled = PETSC_TRUE;
716: C->preallocated = PETSC_TRUE;
717: PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
718: return(0);
719: }
723: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
724: {
725: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ *)C->data;
727: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
728: PetscInt *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
729: PetscInt bs=A->rmap->bs,bs2 = a->bs2,bslog;
730: MatScalar *ba = b->a,*aa,*ap,*dk,*uik;
731: MatScalar *u,*diag,*rtmp,*rtmp_ptr;
732: MatScalar *work;
733: PetscInt *pivots;
736: PetscMalloc(bs2*mbs*sizeof(MatScalar),&rtmp);
737: PetscMemzero(rtmp,bs2*mbs*sizeof(MatScalar));
738: PetscMalloc2(mbs,PetscInt,&il,mbs,PetscInt,&jl);
739: for (i=0; i<mbs; i++) {
740: jl[i] = mbs; il[0] = 0;
741: }
742: PetscMalloc3(bs2,MatScalar,&dk,bs2,MatScalar,&uik,bs,MatScalar,&work);
743: PetscMalloc(bs*sizeof(PetscInt),&pivots);
744:
745: ai = a->i; aj = a->j; aa = a->a;
747: /* flops in while loop */
748: bslog = 2*bs*bs2;
749:
750: /* for each row k */
751: for (k = 0; k<mbs; k++){
753: /*initialize k-th row with elements nonzero in row k of A */
754: jmin = ai[k]; jmax = ai[k+1];
755: ap = aa + jmin*bs2;
756: for (j = jmin; j < jmax; j++){
757: vj = aj[j]; /* block col. index */
758: rtmp_ptr = rtmp + vj*bs2;
759: for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
760: }
762: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
763: PetscMemcpy(dk,rtmp+k*bs2,bs2*sizeof(MatScalar));
764: i = jl[k]; /* first row to be added to k_th row */
766: while (i < k){
767: nexti = jl[i]; /* next row to be added to k_th row */
769: /* compute multiplier */
770: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
772: /* uik = -inv(Di)*U_bar(i,k) */
773: diag = ba + i*bs2;
774: u = ba + ili*bs2;
775: PetscMemzero(uik,bs2*sizeof(MatScalar));
776: PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);
777:
778: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
779: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
780: PetscLogFlops(bslog*2.0);
781:
782: /* update -U(i,k) */
783: PetscMemcpy(ba+ili*bs2,uik,bs2*sizeof(MatScalar));
785: /* add multiple of row i to k-th row ... */
786: jmin = ili + 1; jmax = bi[i+1];
787: if (jmin < jmax){
788: for (j=jmin; j<jmax; j++) {
789: /* rtmp += -U(i,k)^T * U_bar(i,j) */
790: rtmp_ptr = rtmp + bj[j]*bs2;
791: u = ba + j*bs2;
792: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
793: }
794: PetscLogFlops(bslog*(jmax-jmin));
795:
796: /* ... add i to row list for next nonzero entry */
797: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
798: j = bj[jmin];
799: jl[i] = jl[j]; jl[j] = i; /* update jl */
800: }
801: i = nexti;
802: }
804: /* save nonzero entries in k-th row of U ... */
806: /* invert diagonal block */
807: diag = ba+k*bs2;
808: PetscMemcpy(diag,dk,bs2*sizeof(MatScalar));
809: PetscKernel_A_gets_inverse_A(bs,diag,pivots,work);
810:
811: jmin = bi[k]; jmax = bi[k+1];
812: if (jmin < jmax) {
813: for (j=jmin; j<jmax; j++){
814: vj = bj[j]; /* block col. index of U */
815: u = ba + j*bs2;
816: rtmp_ptr = rtmp + vj*bs2;
817: for (k1=0; k1<bs2; k1++){
818: *u++ = *rtmp_ptr;
819: *rtmp_ptr++ = 0.0;
820: }
821: }
822:
823: /* ... add k to row list for first nonzero entry in k-th row */
824: il[k] = jmin;
825: i = bj[jmin];
826: jl[k] = jl[i]; jl[i] = k;
827: }
828: }
830: PetscFree(rtmp);
831: PetscFree2(il,jl);
832: PetscFree3(dk,uik,work);
833: PetscFree(pivots);
835: C->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
836: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
837: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
838: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
839: C->assembled = PETSC_TRUE;
840: C->preallocated = PETSC_TRUE;
841: PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
842: return(0);
843: }
845: /*
846: Numeric U^T*D*U factorization for SBAIJ format. Modified from SNF of YSMP.
847: Version for blocks 2 by 2.
848: */
851: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat C,Mat A,const MatFactorInfo *info)
852: {
853: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ *)C->data;
854: IS perm = b->row;
856: const PetscInt *ai,*aj,*perm_ptr;
857: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
858: PetscInt *a2anew,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
859: MatScalar *ba = b->a,*aa,*ap;
860: MatScalar *u,*diag,*rtmp,*rtmp_ptr,dk[4],uik[4];
861: PetscReal shift = info->shiftamount;
864: /* initialization */
865: /* il and jl record the first nonzero element in each row of the accessing
866: window U(0:k, k:mbs-1).
867: jl: list of rows to be added to uneliminated rows
868: i>= k: jl(i) is the first row to be added to row i
869: i< k: jl(i) is the row following row i in some list of rows
870: jl(i) = mbs indicates the end of a list
871: il(i): points to the first nonzero element in columns k,...,mbs-1 of
872: row i of U */
873: PetscMalloc(4*mbs*sizeof(MatScalar),&rtmp);
874: PetscMemzero(rtmp,4*mbs*sizeof(MatScalar));
875: PetscMalloc2(mbs,PetscInt,&il,mbs,PetscInt,&jl);
876: for (i=0; i<mbs; i++) {
877: jl[i] = mbs; il[0] = 0;
878: }
879: ISGetIndices(perm,&perm_ptr);
881: /* check permutation */
882: if (!a->permute){
883: ai = a->i; aj = a->j; aa = a->a;
884: } else {
885: ai = a->inew; aj = a->jnew;
886: PetscMalloc(4*ai[mbs]*sizeof(MatScalar),&aa);
887: PetscMemcpy(aa,a->a,4*ai[mbs]*sizeof(MatScalar));
888: PetscMalloc(ai[mbs]*sizeof(PetscInt),&a2anew);
889: PetscMemcpy(a2anew,a->a2anew,(ai[mbs])*sizeof(PetscInt));
891: for (i=0; i<mbs; i++){
892: jmin = ai[i]; jmax = ai[i+1];
893: for (j=jmin; j<jmax; j++){
894: while (a2anew[j] != j){
895: k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
896: for (k1=0; k1<4; k1++){
897: dk[k1] = aa[k*4+k1];
898: aa[k*4+k1] = aa[j*4+k1];
899: aa[j*4+k1] = dk[k1];
900: }
901: }
902: /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
903: if (i > aj[j]){
904: /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
905: ap = aa + j*4; /* ptr to the beginning of the block */
906: dk[1] = ap[1]; /* swap ap[1] and ap[2] */
907: ap[1] = ap[2];
908: ap[2] = dk[1];
909: }
910: }
911: }
912: PetscFree(a2anew);
913: }
915: /* for each row k */
916: for (k = 0; k<mbs; k++){
918: /*initialize k-th row with elements nonzero in row perm(k) of A */
919: jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
920: ap = aa + jmin*4;
921: for (j = jmin; j < jmax; j++){
922: vj = perm_ptr[aj[j]]; /* block col. index */
923: rtmp_ptr = rtmp + vj*4;
924: for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
925: }
927: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
928: PetscMemcpy(dk,rtmp+k*4,4*sizeof(MatScalar));
929: i = jl[k]; /* first row to be added to k_th row */
931: while (i < k){
932: nexti = jl[i]; /* next row to be added to k_th row */
934: /* compute multiplier */
935: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
937: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
938: diag = ba + i*4;
939: u = ba + ili*4;
940: uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
941: uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
942: uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
943: uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);
944:
945: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
946: dk[0] += uik[0]*u[0] + uik[1]*u[1];
947: dk[1] += uik[2]*u[0] + uik[3]*u[1];
948: dk[2] += uik[0]*u[2] + uik[1]*u[3];
949: dk[3] += uik[2]*u[2] + uik[3]*u[3];
951: PetscLogFlops(16.0*2.0);
953: /* update -U(i,k): ba[ili] = uik */
954: PetscMemcpy(ba+ili*4,uik,4*sizeof(MatScalar));
956: /* add multiple of row i to k-th row ... */
957: jmin = ili + 1; jmax = bi[i+1];
958: if (jmin < jmax){
959: for (j=jmin; j<jmax; j++) {
960: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
961: rtmp_ptr = rtmp + bj[j]*4;
962: u = ba + j*4;
963: rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
964: rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
965: rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
966: rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
967: }
968: PetscLogFlops(16.0*(jmax-jmin));
969:
970: /* ... add i to row list for next nonzero entry */
971: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
972: j = bj[jmin];
973: jl[i] = jl[j]; jl[j] = i; /* update jl */
974: }
975: i = nexti;
976: }
978: /* save nonzero entries in k-th row of U ... */
980: /* invert diagonal block */
981: diag = ba+k*4;
982: PetscMemcpy(diag,dk,4*sizeof(MatScalar));
983: PetscKernel_A_gets_inverse_A_2(diag,shift);
984:
985: jmin = bi[k]; jmax = bi[k+1];
986: if (jmin < jmax) {
987: for (j=jmin; j<jmax; j++){
988: vj = bj[j]; /* block col. index of U */
989: u = ba + j*4;
990: rtmp_ptr = rtmp + vj*4;
991: for (k1=0; k1<4; k1++){
992: *u++ = *rtmp_ptr;
993: *rtmp_ptr++ = 0.0;
994: }
995: }
996:
997: /* ... add k to row list for first nonzero entry in k-th row */
998: il[k] = jmin;
999: i = bj[jmin];
1000: jl[k] = jl[i]; jl[i] = k;
1001: }
1002: }
1004: PetscFree(rtmp);
1005: PetscFree2(il,jl);
1006: if (a->permute) {
1007: PetscFree(aa);
1008: }
1009: ISRestoreIndices(perm,&perm_ptr);
1010: C->ops->solve = MatSolve_SeqSBAIJ_2_inplace;
1011: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_inplace;
1012: C->assembled = PETSC_TRUE;
1013: C->preallocated = PETSC_TRUE;
1014: PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1015: return(0);
1016: }
1018: /*
1019: Version for when blocks are 2 by 2 Using natural ordering
1020: */
1023: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
1024: {
1025: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ *)C->data;
1027: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
1028: PetscInt *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
1029: MatScalar *ba = b->a,*aa,*ap,dk[8],uik[8];
1030: MatScalar *u,*diag,*rtmp,*rtmp_ptr;
1031: PetscReal shift = info->shiftamount;
1034: /* initialization */
1035: /* il and jl record the first nonzero element in each row of the accessing
1036: window U(0:k, k:mbs-1).
1037: jl: list of rows to be added to uneliminated rows
1038: i>= k: jl(i) is the first row to be added to row i
1039: i< k: jl(i) is the row following row i in some list of rows
1040: jl(i) = mbs indicates the end of a list
1041: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1042: row i of U */
1043: PetscMalloc(4*mbs*sizeof(MatScalar),&rtmp);
1044: PetscMemzero(rtmp,4*mbs*sizeof(MatScalar));
1045: PetscMalloc2(mbs,PetscInt,&il,mbs,PetscInt,&jl);
1046: for (i=0; i<mbs; i++) {
1047: jl[i] = mbs; il[0] = 0;
1048: }
1049: ai = a->i; aj = a->j; aa = a->a;
1051: /* for each row k */
1052: for (k = 0; k<mbs; k++){
1054: /*initialize k-th row with elements nonzero in row k of A */
1055: jmin = ai[k]; jmax = ai[k+1];
1056: ap = aa + jmin*4;
1057: for (j = jmin; j < jmax; j++){
1058: vj = aj[j]; /* block col. index */
1059: rtmp_ptr = rtmp + vj*4;
1060: for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
1061: }
1062:
1063: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
1064: PetscMemcpy(dk,rtmp+k*4,4*sizeof(MatScalar));
1065: i = jl[k]; /* first row to be added to k_th row */
1067: while (i < k){
1068: nexti = jl[i]; /* next row to be added to k_th row */
1070: /* compute multiplier */
1071: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1073: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
1074: diag = ba + i*4;
1075: u = ba + ili*4;
1076: uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
1077: uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
1078: uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
1079: uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);
1080:
1081: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
1082: dk[0] += uik[0]*u[0] + uik[1]*u[1];
1083: dk[1] += uik[2]*u[0] + uik[3]*u[1];
1084: dk[2] += uik[0]*u[2] + uik[1]*u[3];
1085: dk[3] += uik[2]*u[2] + uik[3]*u[3];
1087: PetscLogFlops(16.0*2.0);
1089: /* update -U(i,k): ba[ili] = uik */
1090: PetscMemcpy(ba+ili*4,uik,4*sizeof(MatScalar));
1092: /* add multiple of row i to k-th row ... */
1093: jmin = ili + 1; jmax = bi[i+1];
1094: if (jmin < jmax){
1095: for (j=jmin; j<jmax; j++) {
1096: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
1097: rtmp_ptr = rtmp + bj[j]*4;
1098: u = ba + j*4;
1099: rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
1100: rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
1101: rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
1102: rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
1103: }
1104: PetscLogFlops(16.0*(jmax-jmin));
1106: /* ... add i to row list for next nonzero entry */
1107: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
1108: j = bj[jmin];
1109: jl[i] = jl[j]; jl[j] = i; /* update jl */
1110: }
1111: i = nexti;
1112: }
1114: /* save nonzero entries in k-th row of U ... */
1116: /* invert diagonal block */
1117: diag = ba+k*4;
1118: PetscMemcpy(diag,dk,4*sizeof(MatScalar));
1119: PetscKernel_A_gets_inverse_A_2(diag,shift);
1120:
1121: jmin = bi[k]; jmax = bi[k+1];
1122: if (jmin < jmax) {
1123: for (j=jmin; j<jmax; j++){
1124: vj = bj[j]; /* block col. index of U */
1125: u = ba + j*4;
1126: rtmp_ptr = rtmp + vj*4;
1127: for (k1=0; k1<4; k1++){
1128: *u++ = *rtmp_ptr;
1129: *rtmp_ptr++ = 0.0;
1130: }
1131: }
1132:
1133: /* ... add k to row list for first nonzero entry in k-th row */
1134: il[k] = jmin;
1135: i = bj[jmin];
1136: jl[k] = jl[i]; jl[i] = k;
1137: }
1138: }
1140: PetscFree(rtmp);
1141: PetscFree2(il,jl);
1143: C->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1144: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1145: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1146: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1147: C->assembled = PETSC_TRUE;
1148: C->preallocated = PETSC_TRUE;
1149: PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1150: return(0);
1151: }
1153: /*
1154: Numeric U^T*D*U factorization for SBAIJ format.
1155: Version for blocks are 1 by 1.
1156: */
1159: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
1160: {
1161: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ *)C->data;
1162: IS ip=b->row;
1164: const PetscInt *ai,*aj,*rip;
1165: PetscInt *a2anew,i,j,mbs=a->mbs,*bi=b->i,*bj=b->j,*bcol;
1166: PetscInt k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
1167: MatScalar *rtmp,*ba=b->a,*bval,*aa,dk,uikdi;
1168: PetscReal rs;
1169: FactorShiftCtx sctx;
1172: /* MatPivotSetUp(): initialize shift context sctx */
1173: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1174:
1175: ISGetIndices(ip,&rip);
1176: if (!a->permute){
1177: ai = a->i; aj = a->j; aa = a->a;
1178: } else {
1179: ai = a->inew; aj = a->jnew;
1180: nz = ai[mbs];
1181: PetscMalloc(nz*sizeof(MatScalar),&aa);
1182: a2anew = a->a2anew;
1183: bval = a->a;
1184: for (j=0; j<nz; j++){
1185: aa[a2anew[j]] = *(bval++);
1186: }
1187: }
1188:
1189: /* initialization */
1190: /* il and jl record the first nonzero element in each row of the accessing
1191: window U(0:k, k:mbs-1).
1192: jl: list of rows to be added to uneliminated rows
1193: i>= k: jl(i) is the first row to be added to row i
1194: i< k: jl(i) is the row following row i in some list of rows
1195: jl(i) = mbs indicates the end of a list
1196: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1197: row i of U */
1198: PetscMalloc3(mbs,MatScalar,&rtmp,mbs,PetscInt,&il,mbs,PetscInt,&jl);
1200: do {
1201: sctx.newshift = PETSC_FALSE;
1202: for (i=0; i<mbs; i++) {
1203: rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
1204: }
1205:
1206: for (k = 0; k<mbs; k++){
1207: /*initialize k-th row by the perm[k]-th row of A */
1208: jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
1209: bval = ba + bi[k];
1210: for (j = jmin; j < jmax; j++){
1211: col = rip[aj[j]];
1212: rtmp[col] = aa[j];
1213: *bval++ = 0.0; /* for in-place factorization */
1214: }
1216: /* shift the diagonal of the matrix */
1217: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1219: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1220: dk = rtmp[k];
1221: i = jl[k]; /* first row to be added to k_th row */
1223: while (i < k){
1224: nexti = jl[i]; /* next row to be added to k_th row */
1226: /* compute multiplier, update diag(k) and U(i,k) */
1227: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1228: uikdi = - ba[ili]*ba[bi[i]]; /* diagonal(k) */
1229: dk += uikdi*ba[ili];
1230: ba[ili] = uikdi; /* -U(i,k) */
1232: /* add multiple of row i to k-th row */
1233: jmin = ili + 1; jmax = bi[i+1];
1234: if (jmin < jmax){
1235: for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1236: PetscLogFlops(2.0*(jmax-jmin));
1238: /* update il and jl for row i */
1239: il[i] = jmin;
1240: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1241: }
1242: i = nexti;
1243: }
1245: /* shift the diagonals when zero pivot is detected */
1246: /* compute rs=sum of abs(off-diagonal) */
1247: rs = 0.0;
1248: jmin = bi[k]+1;
1249: nz = bi[k+1] - jmin;
1250: if (nz){
1251: bcol = bj + jmin;
1252: while (nz--){
1253: rs += PetscAbsScalar(rtmp[*bcol]);
1254: bcol++;
1255: }
1256: }
1258: sctx.rs = rs;
1259: sctx.pv = dk;
1260: MatPivotCheck(A,info,&sctx,k);
1261: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1262: dk = sctx.pv;
1263:
1264: /* copy data into U(k,:) */
1265: ba[bi[k]] = 1.0/dk; /* U(k,k) */
1266: jmin = bi[k]+1; jmax = bi[k+1];
1267: if (jmin < jmax) {
1268: for (j=jmin; j<jmax; j++){
1269: col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
1270: }
1271: /* add the k-th row into il and jl */
1272: il[k] = jmin;
1273: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1274: }
1275: }
1276: } while (sctx.newshift);
1277: PetscFree3(rtmp,il,jl);
1278: if (a->permute){PetscFree(aa);}
1280: ISRestoreIndices(ip,&rip);
1281: C->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
1282: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1283: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
1284: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
1285: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
1286: C->assembled = PETSC_TRUE;
1287: C->preallocated = PETSC_TRUE;
1288: PetscLogFlops(C->rmap->N);
1289: if (sctx.nshift){
1290: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1291: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
1292: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1293: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
1294: }
1295: }
1296: return(0);
1297: }
1299: /*
1300: Version for when blocks are 1 by 1 Using natural ordering under new datastructure
1301: Modified from MatCholeskyFactorNumeric_SeqAIJ()
1302: */
1305: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
1306: {
1307: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
1308: Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)B->data;
1310: PetscInt i,j,mbs=A->rmap->n,*bi=b->i,*bj=b->j,*bdiag=b->diag,*bjtmp;
1311: PetscInt *ai=a->i,*aj=a->j,*ajtmp;
1312: PetscInt k,jmin,jmax,*c2r,*il,col,nexti,ili,nz;
1313: MatScalar *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
1314: FactorShiftCtx sctx;
1315: PetscReal rs;
1316: MatScalar d,*v;
1319: PetscMalloc3(mbs,MatScalar,&rtmp,mbs,PetscInt,&il,mbs,PetscInt,&c2r);
1321: /* MatPivotSetUp(): initialize shift context sctx */
1322: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1324: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1325: sctx.shift_top = info->zeropivot;
1326: PetscMemzero(rtmp,mbs*sizeof(MatScalar));
1327: for (i=0; i<mbs; i++) {
1328: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1329: d = (aa)[a->diag[i]];
1330: rtmp[i] += - PetscRealPart(d); /* diagonal entry */
1331: ajtmp = aj + ai[i] + 1; /* exclude diagonal */
1332: v = aa + ai[i] + 1;
1333: nz = ai[i+1] - ai[i] - 1 ;
1334: for (j=0; j<nz; j++){
1335: rtmp[i] += PetscAbsScalar(v[j]);
1336: rtmp[ajtmp[j]] += PetscAbsScalar(v[j]);
1337: }
1338: if (PetscRealPart(rtmp[i]) > sctx.shift_top) sctx.shift_top = PetscRealPart(rtmp[i]);
1339: }
1340: sctx.shift_top *= 1.1;
1341: sctx.nshift_max = 5;
1342: sctx.shift_lo = 0.;
1343: sctx.shift_hi = 1.;
1344: }
1346: /* allocate working arrays
1347: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1348: il: for active k row, il[i] gives the index of the 1st nonzero entry in U[i,k:n-1] in bj and ba arrays
1349: */
1350: do {
1351: sctx.newshift = PETSC_FALSE;
1353: for (i=0; i<mbs; i++) c2r[i] = mbs;
1354: if (mbs) il[0] = 0;
1355:
1356: for (k = 0; k<mbs; k++){
1357: /* zero rtmp */
1358: nz = bi[k+1] - bi[k];
1359: bjtmp = bj + bi[k];
1360: for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
1361:
1362: /* load in initial unfactored row */
1363: bval = ba + bi[k];
1364: jmin = ai[k]; jmax = ai[k+1];
1365: for (j = jmin; j < jmax; j++){
1366: col = aj[j];
1367: rtmp[col] = aa[j];
1368: *bval++ = 0.0; /* for in-place factorization */
1369: }
1370: /* shift the diagonal of the matrix: ZeropivotApply() */
1371: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1372:
1373: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1374: dk = rtmp[k];
1375: i = c2r[k]; /* first row to be added to k_th row */
1377: while (i < k){
1378: nexti = c2r[i]; /* next row to be added to k_th row */
1379:
1380: /* compute multiplier, update diag(k) and U(i,k) */
1381: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1382: uikdi = - ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
1383: dk += uikdi*ba[ili]; /* update diag[k] */
1384: ba[ili] = uikdi; /* -U(i,k) */
1386: /* add multiple of row i to k-th row */
1387: jmin = ili + 1; jmax = bi[i+1];
1388: if (jmin < jmax){
1389: for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1390: /* update il and c2r for row i */
1391: il[i] = jmin;
1392: j = bj[jmin]; c2r[i] = c2r[j]; c2r[j] = i;
1393: }
1394: i = nexti;
1395: }
1397: /* copy data into U(k,:) */
1398: rs = 0.0;
1399: jmin = bi[k]; jmax = bi[k+1]-1;
1400: if (jmin < jmax) {
1401: for (j=jmin; j<jmax; j++){
1402: col = bj[j]; ba[j] = rtmp[col]; rs += PetscAbsScalar(ba[j]);
1403: }
1404: /* add the k-th row into il and c2r */
1405: il[k] = jmin;
1406: i = bj[jmin]; c2r[k] = c2r[i]; c2r[i] = k;
1407: }
1409: sctx.rs = rs;
1410: sctx.pv = dk;
1411: MatPivotCheck(A,info,&sctx,k);
1412: if(sctx.newshift) break;
1413: dk = sctx.pv;
1414:
1415: ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
1416: }
1417: } while (sctx.newshift);
1418:
1419: PetscFree3(rtmp,il,c2r);
1420:
1421: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1422: B->ops->solves = MatSolves_SeqSBAIJ_1;
1423: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1424: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1425: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1427: B->assembled = PETSC_TRUE;
1428: B->preallocated = PETSC_TRUE;
1429: PetscLogFlops(B->rmap->n);
1431: /* MatPivotView() */
1432: if (sctx.nshift){
1433: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1434: PetscInfo4(A,"number of shift_pd tries %D, shift_amount %G, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,sctx.shift_amount,sctx.shift_fraction,sctx.shift_top);
1435: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1436: PetscInfo2(A,"number of shift_nz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
1437: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS){
1438: PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %G\n",sctx.nshift,info->shiftamount);
1439: }
1440: }
1441: return(0);
1442: }
1446: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
1447: {
1448: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ *)C->data;
1450: PetscInt i,j,mbs = a->mbs;
1451: PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
1452: PetscInt k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
1453: MatScalar *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
1454: PetscReal rs;
1455: FactorShiftCtx sctx;
1458: /* MatPivotSetUp(): initialize shift context sctx */
1459: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1461: /* initialization */
1462: /* il and jl record the first nonzero element in each row of the accessing
1463: window U(0:k, k:mbs-1).
1464: jl: list of rows to be added to uneliminated rows
1465: i>= k: jl(i) is the first row to be added to row i
1466: i< k: jl(i) is the row following row i in some list of rows
1467: jl(i) = mbs indicates the end of a list
1468: il(i): points to the first nonzero element in U(i,k:mbs-1)
1469: */
1470: PetscMalloc(mbs*sizeof(MatScalar),&rtmp);
1471: PetscMalloc2(mbs,PetscInt,&il,mbs,PetscInt,&jl);
1473: do {
1474: sctx.newshift = PETSC_FALSE;
1475: for (i=0; i<mbs; i++) {
1476: rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
1477: }
1479: for (k = 0; k<mbs; k++){
1480: /*initialize k-th row with elements nonzero in row perm(k) of A */
1481: nz = ai[k+1] - ai[k];
1482: acol = aj + ai[k];
1483: aval = aa + ai[k];
1484: bval = ba + bi[k];
1485: while (nz -- ){
1486: rtmp[*acol++] = *aval++;
1487: *bval++ = 0.0; /* for in-place factorization */
1488: }
1490: /* shift the diagonal of the matrix */
1491: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1492:
1493: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1494: dk = rtmp[k];
1495: i = jl[k]; /* first row to be added to k_th row */
1497: while (i < k){
1498: nexti = jl[i]; /* next row to be added to k_th row */
1499: /* compute multiplier, update D(k) and U(i,k) */
1500: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1501: uikdi = - ba[ili]*ba[bi[i]];
1502: dk += uikdi*ba[ili];
1503: ba[ili] = uikdi; /* -U(i,k) */
1505: /* add multiple of row i to k-th row ... */
1506: jmin = ili + 1;
1507: nz = bi[i+1] - jmin;
1508: if (nz > 0){
1509: bcol = bj + jmin;
1510: bval = ba + jmin;
1511: PetscLogFlops(2.0*nz);
1512: while (nz --) rtmp[*bcol++] += uikdi*(*bval++);
1513:
1514: /* update il and jl for i-th row */
1515: il[i] = jmin;
1516: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1517: }
1518: i = nexti;
1519: }
1521: /* shift the diagonals when zero pivot is detected */
1522: /* compute rs=sum of abs(off-diagonal) */
1523: rs = 0.0;
1524: jmin = bi[k]+1;
1525: nz = bi[k+1] - jmin;
1526: if (nz){
1527: bcol = bj + jmin;
1528: while (nz--){
1529: rs += PetscAbsScalar(rtmp[*bcol]);
1530: bcol++;
1531: }
1532: }
1534: sctx.rs = rs;
1535: sctx.pv = dk;
1536: MatPivotCheck(A,info,&sctx,k);
1537: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1538: dk = sctx.pv;
1539:
1540: /* copy data into U(k,:) */
1541: ba[bi[k]] = 1.0/dk;
1542: jmin = bi[k]+1;
1543: nz = bi[k+1] - jmin;
1544: if (nz){
1545: bcol = bj + jmin;
1546: bval = ba + jmin;
1547: while (nz--){
1548: *bval++ = rtmp[*bcol];
1549: rtmp[*bcol++] = 0.0;
1550: }
1551: /* add k-th row into il and jl */
1552: il[k] = jmin;
1553: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1554: }
1555: } /* end of for (k = 0; k<mbs; k++) */
1556: } while (sctx.newshift);
1557: PetscFree(rtmp);
1558: PetscFree2(il,jl);
1559:
1560: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1561: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1562: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1563: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1564: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1566: C->assembled = PETSC_TRUE;
1567: C->preallocated = PETSC_TRUE;
1568: PetscLogFlops(C->rmap->N);
1569: if (sctx.nshift){
1570: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1571: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
1572: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1573: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
1574: }
1575: }
1576: return(0);
1577: }
1581: PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A,IS perm,const MatFactorInfo *info)
1582: {
1584: Mat C;
1587: MatGetFactor(A,"petsc",MAT_FACTOR_CHOLESKY,&C);
1588: MatCholeskyFactorSymbolic(C,A,perm,info);
1589: MatCholeskyFactorNumeric(C,A,info);
1590: A->ops->solve = C->ops->solve;
1591: A->ops->solvetranspose = C->ops->solvetranspose;
1592: MatHeaderMerge(A,C);
1593: return(0);
1594: }