2: #include <../src/mat/utils/freespace.h>

  6: PetscErrorCode PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList *list)
  7: {
  8:   PetscFreeSpaceList a;
  9:   PetscErrorCode     ierr;

 12:   PetscMalloc(sizeof(struct _Space),&a);
 13:   PetscMalloc(n*sizeof(PetscInt),&(a->array_head));
 14:   a->array            = a->array_head;
 15:   a->local_remaining  = n;
 16:   a->local_used       = 0;
 17:   a->total_array_size = 0;
 18:   a->more_space       = PETSC_NULL;

 20:   if (*list) {
 21:     (*list)->more_space = a;
 22:     a->total_array_size = (*list)->total_array_size;
 23:   }

 25:   a->total_array_size += n;
 26:   *list               =  a;
 27:   return(0);
 28: }

 32: PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head,PetscInt *space)
 33: {
 34:   PetscFreeSpaceList a;
 35:   PetscErrorCode     ierr;

 38:   while ((*head)) {
 39:     a     =  (*head)->more_space;
 40:      PetscMemcpy(space,(*head)->array_head,((*head)->local_used)*sizeof(PetscInt));
 41:     space += (*head)->local_used;
 42:      PetscFree((*head)->array_head);
 43:      PetscFree(*head);
 44:     *head =  a;
 45:   }
 46:   return(0);
 47: }

 49: /*
 50:   PetscFreeSpaceContiguous_LU -
 51:     Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array 
 52:   that enables an efficient matrix triangular solve.

 54:    Input Parameters:
 55: +  head - linked list of column indices obtained from matrix symbolic ILU or LU factorization
 56: .  space - an allocated int array with length nnz of factored matrix. 
 57: .  n - order of the matrix
 58: .  bi - row pointer of factored matrix L with length n+1.
 59: -  bdiag - int array of length n+1. bdiag[i] points to diagonal of U(i,:), and bdiag[n] points to entry of U(n-1,0)-1.

 61:    Output Parameter:
 62: .  space - column indices are copied into this int array with contiguous layout of L and U

 64:    See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
 65: */
 68: PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag)
 69: {
 70:   PetscFreeSpaceList a;
 71:   PetscErrorCode     ierr;
 72:   PetscInt           row,nnz,*bj,*array,total,bi_temp;
 73:   PetscInt           nnzL,nnzU;

 76:   bi_temp = bi[n];
 77:   row       = 0;
 78:   total     = 0;
 79:   nnzL  = bdiag[0];
 80:   while ((*head)) {
 81:     total += (*head)->local_used;
 82:     array  = (*head)->array_head;
 83: 
 84:     while (row < n) {
 85:       if (bi[row+1] > total) break;
 86:       /* copy array entries into bj for this row */
 87:       nnz  = bi[row+1] - bi[row];
 88:       /* set bi[row] for new datastruct */
 89:       if (row == 0 ){
 90:         bi[row] = 0;
 91:       } else {
 92:         bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */
 93:       }

 95:       /* L part */
 96:       nnzL = bdiag[row];
 97:       bj   = space+bi[row];
 98:       PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));
 99: 
100:       /* diagonal entry */
101:       bdiag[row] = bi_temp - 1;
102:       space[bdiag[row]] = row;

104:       /* U part */
105:       nnzU        = nnz - nnzL;
106:       bi_temp = bi_temp - nnzU;
107:       nnzU --;      /* exclude diagonal */
108:       bj = space + bi_temp;
109:       PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));
110:       array += nnz;
111:       row++;
112:     }

114:     a     = (*head)->more_space;
115:     PetscFree((*head)->array_head);
116:     PetscFree(*head);
117:     *head = a;
118:   }
119:   if (n) {
120:     bi[n] = bi[n-1] + nnzL;
121:     bdiag[n] = bdiag[n-1]-1;
122:   }
123:   return(0);
124: }

126: /*
127:   PetscFreeSpaceContiguous_Cholesky -
128:     Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array 
129:   that enables an efficient matrix triangular solve.

131:    Input Parameters:
132: +  head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization
133: .  space - an allocated int array with length nnz of factored matrix. 
134: .  n - order of the matrix
135: .  ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix.
136: -  udiag - int array of length n.

138:    Output Parameter:
139: +  space - column indices are copied into this int array with contiguous layout of U, with diagonal located as the last entry in each row
140: -  udiag - indices of diagonal entries

142:    See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
143: */

147: PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag)
148: {
149:   PetscFreeSpaceList a;
150:   PetscErrorCode     ierr;
151:   PetscInt           row,nnz,*uj,*array,total;

154:   row   = 0;
155:   total = 0;
156:   while (*head) {
157:     total += (*head)->local_used;
158:     array  = (*head)->array_head;
159: 
160:     while (row < n){
161:       if (ui[row+1] > total) break;
162:       udiag[row] = ui[row+1] - 1;     /* points to the last entry of U(row,:) */
163:       nnz  = ui[row+1] - ui[row] - 1; /* exclude diagonal */
164:       uj   = space + ui[row];
165:       PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));
166:       uj[nnz] = array[0]; /* diagonal */
167:       array += nnz + 1;
168:       row++;
169:     }

171:     a     = (*head)->more_space;
172:     PetscFree((*head)->array_head);
173:     PetscFree(*head);
174:     *head = a;
175:   }
176:   return(0);
177: }

181: PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
182: {
183:   PetscFreeSpaceList a;
184:   PetscErrorCode     ierr;

187:   while ((head)) {
188:     a    = (head)->more_space;
189:     PetscFree((head)->array_head);
190:     PetscFree(head);
191:     head = a;
192:   }
193:   return(0);
194: }