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