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 15 a->array = a->array_head; 16 a->local_remaining = n; 17 a->local_used = 0; 18 a->total_array_size = 0; 19 a->more_space = NULL; 20 21 if (*list) { 22 (*list)->more_space = a; 23 a->total_array_size = (*list)->total_array_size; 24 } 25 26 a->total_array_size += n; 27 *list = a; 28 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "PetscFreeSpaceContiguous" 33 PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head,PetscInt *space) 34 { 35 PetscFreeSpaceList a; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 while ((*head)) { 40 a = (*head)->more_space; 41 ierr = PetscMemcpy(space,(*head)->array_head,((*head)->local_used)*sizeof(PetscInt));CHKERRQ(ierr); 42 space += (*head)->local_used; 43 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 44 ierr = PetscFree(*head);CHKERRQ(ierr); 45 *head = a; 46 } 47 PetscFunctionReturn(0); 48 } 49 50 /* 51 PetscFreeSpaceContiguous_LU - 52 Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array 53 that enables an efficient matrix triangular solve. 54 55 Input Parameters: 56 + head - linked list of column indices obtained from matrix symbolic ILU or LU factorization 57 . space - an allocated int array with length nnz of factored matrix. 58 . n - order of the matrix 59 . bi - row pointer of factored matrix L with length n+1. 60 - 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 62 Output Parameter: 63 . space - column indices are copied into this int array with contiguous layout of L and U 64 65 See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U 66 */ 67 #undef __FUNCT__ 68 #define __FUNCT__ "PetscFreeSpaceContiguous_LU" 69 PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag) 70 { 71 PetscFreeSpaceList a; 72 PetscErrorCode ierr; 73 PetscInt row,nnz,*bj,*array,total,bi_temp; 74 PetscInt nnzL,nnzU; 75 76 PetscFunctionBegin; 77 bi_temp = bi[n]; 78 row = 0; 79 total = 0; 80 nnzL = bdiag[0]; 81 while ((*head)) { 82 total += (*head)->local_used; 83 array = (*head)->array_head; 84 85 while (row < n) { 86 if (bi[row+1] > total) break; 87 /* copy array entries into bj for this row */ 88 nnz = bi[row+1] - bi[row]; 89 /* set bi[row] for new datastruct */ 90 if (row == 0) { 91 bi[row] = 0; 92 } else { 93 bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */ 94 } 95 96 /* L part */ 97 nnzL = bdiag[row]; 98 bj = space+bi[row]; 99 ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr); 100 101 /* diagonal entry */ 102 bdiag[row] = bi_temp - 1; 103 space[bdiag[row]] = row; 104 105 /* U part */ 106 nnzU = nnz - nnzL; 107 bi_temp = bi_temp - nnzU; 108 nnzU--; /* exclude diagonal */ 109 bj = space + bi_temp; 110 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 111 array += nnz; 112 row++; 113 } 114 115 a = (*head)->more_space; 116 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 117 ierr = PetscFree(*head);CHKERRQ(ierr); 118 *head = a; 119 } 120 if (n) { 121 bi[n] = bi[n-1] + nnzL; 122 bdiag[n] = bdiag[n-1]-1; 123 } 124 PetscFunctionReturn(0); 125 } 126 127 /* 128 PetscFreeSpaceContiguous_Cholesky - 129 Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array 130 that enables an efficient matrix triangular solve. 131 132 Input Parameters: 133 + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization 134 . space - an allocated int array with length nnz of factored matrix. 135 . n - order of the matrix 136 . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix. 137 - udiag - int array of length n. 138 139 Output Parameter: 140 + space - column indices are copied into this int array with contiguous layout of U, with diagonal located as the last entry in each row 141 - udiag - indices of diagonal entries 142 143 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description. 144 */ 145 146 #undef __FUNCT__ 147 #define __FUNCT__ "PetscFreeSpaceContiguous_Cholesky" 148 PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag) 149 { 150 PetscFreeSpaceList a; 151 PetscErrorCode ierr; 152 PetscInt row,nnz,*uj,*array,total; 153 154 PetscFunctionBegin; 155 row = 0; 156 total = 0; 157 while (*head) { 158 total += (*head)->local_used; 159 array = (*head)->array_head; 160 161 while (row < n) { 162 if (ui[row+1] > total) break; 163 udiag[row] = ui[row+1] - 1; /* points to the last entry of U(row,:) */ 164 nnz = ui[row+1] - ui[row] - 1; /* exclude diagonal */ 165 uj = space + ui[row]; 166 ierr = PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));CHKERRQ(ierr); 167 uj[nnz] = array[0]; /* diagonal */ 168 array += nnz + 1; 169 row++; 170 } 171 172 a = (*head)->more_space; 173 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 174 ierr = PetscFree(*head);CHKERRQ(ierr); 175 *head = a; 176 } 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "PetscFreeSpaceDestroy" 182 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head) 183 { 184 PetscFreeSpaceList a; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 while ((head)) { 189 a = (head)->more_space; 190 ierr = PetscFree((head)->array_head);CHKERRQ(ierr); 191 ierr = PetscFree(head);CHKERRQ(ierr); 192 head = a; 193 } 194 PetscFunctionReturn(0); 195 } 196