xref: /petsc/src/mat/utils/freespace.c (revision 7319c6545cdf078a8212d19c3cdefe116c1a5fdd)
1 #define PETSCMAT_DLL
2 
3 #include "../src/mat/utils/freespace.h"
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "PetscFreeSpaceGet"
7 PetscErrorCode PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList *list)
8 {
9   PetscFreeSpaceList a;
10   PetscErrorCode     ierr;
11 
12   PetscFunctionBegin;
13   ierr = PetscMalloc(sizeof(struct _Space),&a);CHKERRQ(ierr);
14   ierr = PetscMalloc(n*sizeof(PetscInt),&(a->array_head));CHKERRQ(ierr);
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)!=NULL) {
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   Copy a linket list obtained from matrix symbolic ilu or lu factorization into a contiguous array that enables
52   an efficient matrix 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 arry with length nnz of factored matrix.
57 .  n - order of the matrix
58 .  bi - row pointer of factored matrix with length 2n+2. First n+1 entries are set based on the traditional layout of L and U matrices
59 .
60 -  bdiag - int array holding the number of nonzeros in each row of L matrix, excluding diagonals.
61 
62    Output Parameter:
63 .  space - column indices are copied into this int array with contiguous layout of L and U
64    See MatILUFactorSymbolic_SeqAIJ_ilu0_newdatastruct() for detailed description.
65 */
66 #undef __FUNCT__
67 #define __FUNCT__ "PetscFreeSpaceContiguous_newdatastruct"
68 PetscErrorCode PetscFreeSpaceContiguous_newdatastruct(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag)
69 {
70   PetscFreeSpaceList a;
71   PetscErrorCode     ierr;
72   PetscInt           row,nnz,*bj,*array,total;
73   PetscInt           nnzL,nnzU;
74 
75   PetscFunctionBegin;
76   bi[2*n+1] = bi[n];
77   row       = 0;
78   total     = 0;
79   nnzL  = bdiag[0];
80   while ((*head)!=NULL) {
81     total += (*head)->local_used;
82     array  = (*head)->array_head;
83 
84     while (bi[row+1] <= total && row < n){
85       /* copy array entries into bj for this row */
86       nnz  = bi[row+1] - bi[row];
87       /* set bi[row] for new datastruct */
88       if (row == 0 ){
89         bi[row] = 0;
90       } else {
91         bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */
92       }
93 
94       /* L part */
95       nnzL = bdiag[row];
96       bj   = space+bi[row];
97       ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr);
98 
99       /* diagonal entry */
100       bdiag[row]        = bi[2*n-row+1]-1;
101       space[bdiag[row]] = row;
102 
103       /* U part */
104       nnzU        = nnz - nnzL;
105       bi[2*n-row] = bi[2*n-row+1] - nnzU;
106       nnzU --;      /* exclude diagonal */
107       bj   = space + bi[2*n-(row)];
108       ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr);
109 
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   bi[n] = bi[n-1] + nnzL;
120   if (bi[n] != bi[n+1]) SETERRQ2(1,"bi[n] %d != bi[n+1] %d",bi[n],bi[n+1]);
121   PetscFunctionReturn(0);
122 }
123 
124 #undef __FUNCT__
125 #define __FUNCT__ "PetscFreeSpaceDestroy"
126 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
127 {
128   PetscFreeSpaceList a;
129   PetscErrorCode     ierr;
130 
131   PetscFunctionBegin;
132   while ((head)!=NULL) {
133     a    = (head)->more_space;
134     ierr = PetscFree((head)->array_head);CHKERRQ(ierr);
135     ierr = PetscFree(head);CHKERRQ(ierr);
136     head = a;
137   }
138   PetscFunctionReturn(0);
139 }
140