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