1 2 #include <../src/mat/utils/freespace.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "PetscFreeSpaceGet" 6 PetscErrorCode PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList *list) 7 { 8 PetscFreeSpaceList a; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = PetscMalloc(sizeof(struct _Space),&a);CHKERRQ(ierr); 13 ierr = PetscMalloc(n*sizeof(PetscInt),&(a->array_head));CHKERRQ(ierr); 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; 19 20 if (*list) { 21 (*list)->more_space = a; 22 a->total_array_size = (*list)->total_array_size; 23 } 24 25 a->total_array_size += n; 26 *list = a; 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "PetscFreeSpaceContiguous" 32 PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head,PetscInt *space) 33 { 34 PetscFreeSpaceList a; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 while ((*head)) { 39 a = (*head)->more_space; 40 ierr = PetscMemcpy(space,(*head)->array_head,((*head)->local_used)*sizeof(PetscInt));CHKERRQ(ierr); 41 space += (*head)->local_used; 42 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 43 ierr = PetscFree(*head);CHKERRQ(ierr); 44 *head = a; 45 } 46 PetscFunctionReturn(0); 47 } 48 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. 53 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. 60 61 Output Parameter: 62 . space - column indices are copied into this int array with contiguous layout of L and U 63 64 See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U 65 */ 66 #undef __FUNCT__ 67 #define __FUNCT__ "PetscFreeSpaceContiguous_LU" 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; 74 75 PetscFunctionBegin; 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 } 94 95 /* L part */ 96 nnzL = bdiag[row]; 97 bj = space+bi[row]; 98 ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr); 99 100 /* diagonal entry */ 101 bdiag[row] = bi_temp - 1; 102 space[bdiag[row]] = row; 103 104 /* U part */ 105 nnzU = nnz - nnzL; 106 bi_temp = bi_temp - nnzU; 107 nnzU --; /* exclude diagonal */ 108 bj = space + bi_temp; 109 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 110 array += nnz; 111 row++; 112 } 113 114 a = (*head)->more_space; 115 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 116 ierr = PetscFree(*head);CHKERRQ(ierr); 117 *head = a; 118 } 119 if (n) { 120 bi[n] = bi[n-1] + nnzL; 121 bdiag[n] = bdiag[n-1]-1; 122 } 123 PetscFunctionReturn(0); 124 } 125 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. 130 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. 137 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 141 142 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description. 143 */ 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "PetscFreeSpaceContiguous_Cholesky" 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; 152 153 PetscFunctionBegin; 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 ierr = PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));CHKERRQ(ierr); 166 uj[nnz] = array[0]; /* diagonal */ 167 array += nnz + 1; 168 row++; 169 } 170 171 a = (*head)->more_space; 172 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 173 ierr = PetscFree(*head);CHKERRQ(ierr); 174 *head = a; 175 } 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "PetscFreeSpaceDestroy" 181 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head) 182 { 183 PetscFreeSpaceList a; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 while ((head)) { 188 a = (head)->more_space; 189 ierr = PetscFree((head)->array_head);CHKERRQ(ierr); 190 ierr = PetscFree(head);CHKERRQ(ierr); 191 head = a; 192 } 193 PetscFunctionReturn(0); 194 } 195