1 /*
2  *
3  *  This file is part of MUMPS 5.1.2, released
4  *  on Mon Oct  2 07:37:01 UTC 2017
5  *
6  *
7  *  Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
8  *  University of Bordeaux.
9  *
10  *  This version of MUMPS is provided to you free of charge. It is
11  *  released under the CeCILL-C license:
12  *  http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
13  *
14  */
15 /*
16  * This file contains interfaces to external ordering packages.
17  * At the moment, PORD (J. Schulze) and SCOTCH are interfaced.
18  */
19 #include "mumps_pord.h"
MUMPS_PORD_INTSIZE(MUMPS_INT * pord_intsize)20 void MUMPS_CALL MUMPS_PORD_INTSIZE(MUMPS_INT *pord_intsize)
21 {
22 #if defined(pord)
23 #   if defined(PORD_INTSIZE64) || defined(INTSIZE64)
24     *pord_intsize=64;
25 #   else
26     *pord_intsize=32;
27 #   endif
28 #else
29     *pord_intsize=-99999;
30 #endif
31 }
32 #if defined(pord)
33 /* Interface to PORD */
34 #if defined(INTSIZE64) || defined(PORD_INTSIZE64)
35 void MUMPS_CALL
MUMPS_PORDF(MUMPS_INT8 * nvtx,MUMPS_INT8 * nedges,MUMPS_INT8 * xadj,MUMPS_INT8 * adjncy,MUMPS_INT8 * nv,MUMPS_INT * ncmpa)36 MUMPS_PORDF( MUMPS_INT8 *nvtx, MUMPS_INT8 *nedges,
37              MUMPS_INT8 *xadj, MUMPS_INT8 *adjncy,
38              MUMPS_INT8 *nv, MUMPS_INT *ncmpa )
39 #else
40 void MUMPS_CALL
41 MUMPS_PORDF( MUMPS_INT *nvtx, MUMPS_INT *nedges,
42              MUMPS_INT *xadj, MUMPS_INT *adjncy,
43              MUMPS_INT *nv, MUMPS_INT *ncmpa )
44 #endif
45 {
46     *ncmpa = mumps_pord( *nvtx, *nedges, xadj, adjncy, nv );
47 }
48 /* Interface to PORD with weighted graph */
49 #if defined(INTSIZE64) || defined(PORD_INTSIZE64)
50 void MUMPS_CALL
MUMPS_PORDF_WND(MUMPS_INT8 * nvtx,MUMPS_INT8 * nedges,MUMPS_INT8 * xadj,MUMPS_INT8 * adjncy,MUMPS_INT8 * nv,MUMPS_INT * ncmpa,MUMPS_INT8 * totw)51 MUMPS_PORDF_WND( MUMPS_INT8 *nvtx, MUMPS_INT8 *nedges,
52                  MUMPS_INT8 *xadj, MUMPS_INT8 *adjncy,
53                  MUMPS_INT8 *nv, MUMPS_INT *ncmpa, MUMPS_INT8 *totw )
54 #else
55 void MUMPS_CALL
56 MUMPS_PORDF_WND( MUMPS_INT *nvtx, MUMPS_INT *nedges,
57                  MUMPS_INT *xadj, MUMPS_INT *adjncy,
58                  MUMPS_INT *nv, MUMPS_INT *ncmpa, MUMPS_INT *totw )
59 #endif
60 {
61     *ncmpa = mumps_pord_wnd( *nvtx, *nedges, xadj, adjncy, nv, totw );
62 }
63 /************************************************************
64  mumps_pord is used in ana_aux.F
65         permutation and inverse permutation not set on output,
66         but may be printed in default file: "perm_pord" and "iperm_pord",
67         if associated part uncommneted.
68         But, if uncommetnted a bug occurs in psl_ma41_analysi.F
69 ******************************************************************/
70 /*********************************************************/
mumps_pord(PORD_INT nvtx,PORD_INT nedges,PORD_INT * xadj_pe,PORD_INT * adjncy,PORD_INT * nv)71 MUMPS_INT mumps_pord
72 (
73    PORD_INT nvtx,
74    PORD_INT nedges,      /* NZ-like */
75    PORD_INT *xadj_pe,    /* NZ-like */
76    PORD_INT *adjncy,
77    PORD_INT *nv
78 )
79 {
80 /**********************************
81 Arguments:
82 input:
83 -----
84 - nvtx          : dimension of the Problem (N)
85 - nedges        : number of entries (NZ)
86 - adjncy        : non-zeros entries (IW input)
87 input/output:
88 -------------
89 - xadj_pe       : in: pointer through beginning of column non-zeros entries
90                   out: "father array" (PE)
91 ouput:
92 ------
93 - nv            : "nfront array" (NV)
94 *************************************/
95   graph_t    *G;
96   elimtree_t *T;
97   timings_t  cpus[12];
98   options_t  options[] = { SPACE_ORDTYPE, SPACE_NODE_SELECTION1,
99                     SPACE_NODE_SELECTION2, SPACE_NODE_SELECTION3,
100                     SPACE_DOMAIN_SIZE, 0 };
101   PORD_INT *ncolfactor, *ncolupdate, *parent, *vtx2front;
102   PORD_INT *first, *link, nfronts, J, K, u, vertex, vertex_root, count;
103  /* Explicit shifting of indices to be optimized */
104   for (u = nvtx; u >= 0; u--)
105    {
106      xadj_pe[u] = xadj_pe[u] - 1;
107    }
108    for (K = nedges-1; K >= 0; K--)
109    {
110       adjncy[K] = adjncy[K] - 1;
111    }
112    /* initialization of the graph */
113    mymalloc(G, 1, graph_t);
114    G->xadj   = xadj_pe;
115    G->adjncy = adjncy;
116    G->nvtx = nvtx;
117    G->nedges = nedges;
118    /* FIXME: G->vwght and G->tocwght accessed if G->type==UNWEIGHTED? */
119    mymalloc(G->vwght, nvtx, PORD_INT);
120    G->type = UNWEIGHTED;
121    G->totvwght = nvtx;
122    for (u = 0; u < nvtx; u++)
123      G->vwght[u] = 1;
124    /* main function of the Ordering */
125    T = SPACE_ordering(G, options, cpus);
126    nfronts = T->nfronts;
127    ncolfactor = T->ncolfactor;
128    ncolupdate = T->ncolupdate;
129    parent = T->parent;
130  /*    firstchild = T->firstchild; */
131    vtx2front = T->vtx2front;
132     /* -----------------------------------------------------------
133      store the vertices/columns of a front in a bucket structure
134      ----------------------------------------------------------- */
135    mymalloc(first, nfronts, PORD_INT);
136    mymalloc(link, nvtx, PORD_INT);
137    for (J = 0; J < nfronts; J++)
138       first[J] = -1;
139    for (u = nvtx-1; u >= 0; u--)
140       {
141         J = vtx2front[u];
142         link[u] = first[J];
143         first[J] = u;
144       }
145   /* -----------------------------------------------------------
146      fill the two arrays corresponding to the MUMPS tree structure
147      ----------------------------------------------------------- */
148    count = 0;
149    for (K = firstPostorder(T); K != -1; K = nextPostorder(T, K))
150      {
151        vertex_root = first[K];
152        if (vertex_root == -1)
153          {
154            /* Should never happen */
155 #          if defined(PORD_INTSIZE64) || defined(INTSIZE64)
156            printf(" Internal error in mumps_pord, %ld\n",K);
157 #          else
158            printf(" Internal error in mumps_pord, %d\n",K);
159 #          endif
160            exit(-1);
161          }
162        /* for the principal column of the supervariable */
163        if (parent[K] == -1)
164          xadj_pe[vertex_root] = 0; /* root of the tree */
165        else
166          xadj_pe[vertex_root] = - (first[parent[K]]+1);
167        nv[vertex_root] = ncolfactor[K] + ncolupdate[K];
168        count++;
169        for (vertex = link[vertex_root]; vertex != -1; vertex = link[vertex])
170         /* for the secondary columns of the supervariable */
171          {
172            xadj_pe[vertex] = - (vertex_root+1);
173            nv[vertex] = 0;
174           count++;
175          }
176      }
177   /* ----------------------
178      free memory and return
179      ---------------------- */
180   free(first); free(link);
181   free(G->vwght);
182   free(G);
183   freeElimTree(T);
184   return (0);
185 }
186 /*********************************************************/
mumps_pord_wnd(PORD_INT nvtx,PORD_INT nedges,PORD_INT * xadj_pe,PORD_INT * adjncy,PORD_INT * nv,PORD_INT * totw)187 MUMPS_INT mumps_pord_wnd
188 (
189         PORD_INT nvtx,
190         PORD_INT nedges,
191         PORD_INT *xadj_pe,
192         PORD_INT *adjncy,
193         PORD_INT *nv,
194         PORD_INT *totw
195 )
196 {
197 /**********************************
198 Arguments:
199 input:
200 -----
201 - nvtx          : dimension of the Problem (N)
202 - nedges        : number of entries (NZ)
203 - adjncy        : non-zeros entries (IW input)
204 - totw          : sum of the weigth of the vertices
205 input/output:
206 -------------
207 - xadj_pe       : in: pointer through beginning of column non-zeros entries
208                   out: "father array" (PE)
209 - nv            : in: weight of the vertices
210                   out: "nfront array" (NV)
211 *************************************/
212         graph_t    *G;
213         elimtree_t *T;
214         timings_t  cpus[12];
215         options_t  options[] = { SPACE_ORDTYPE, SPACE_NODE_SELECTION1,
216                     SPACE_NODE_SELECTION2, SPACE_NODE_SELECTION3,
217                     SPACE_DOMAIN_SIZE, 0 };
218         PORD_INT *ncolfactor, *ncolupdate, *parent, *vtx2front;
219         PORD_INT *first, *link, nfronts, J, K, u, vertex, vertex_root, count;
220  /* Explicit shifting of indices to be optimized */
221         for (u = nvtx; u >= 0; u--)
222         {
223           xadj_pe[u] = xadj_pe[u] - 1;
224         }
225         for (K = nedges-1; K >= 0; K--)
226         {
227           adjncy[K] = adjncy[K] - 1;
228         }
229         /* initialization of the graph */
230         mymalloc(G, 1, graph_t);
231         G->xadj  = xadj_pe;
232         G->adjncy= adjncy;
233         G->nvtx = nvtx;
234         G->nedges = nedges;
235         G->type = WEIGHTED;
236         G->totvwght = (*totw);
237         /* FIXME: avoid allocation and do: G->vwght=nv; instead? */
238         mymalloc(G->vwght, nvtx, PORD_INT);
239         for (u = 0; u < nvtx; u++)
240           G->vwght[u] = nv[u];
241         /* main function of the Ordering */
242         T = SPACE_ordering(G, options, cpus);
243         nfronts = T->nfronts;
244         ncolfactor = T->ncolfactor;
245         ncolupdate = T->ncolupdate;
246         parent = T->parent;
247   /*    firstchild = T->firstchild; */
248         vtx2front = T->vtx2front;
249     /* -----------------------------------------------------------
250      store the vertices/columns of a front in a bucket structure
251      ----------------------------------------------------------- */
252         mymalloc(first, nfronts, PORD_INT);
253         mymalloc(link, nvtx, PORD_INT);
254         for (J = 0; J < nfronts; J++)
255           first[J] = -1;
256         for (u = nvtx-1; u >= 0; u--)
257         {
258           J = vtx2front[u];
259           link[u] = first[J];
260           first[J] = u;
261         }
262   /* -----------------------------------------------------------
263      fill the two arrays corresponding to the MUMPS tree structure
264      ----------------------------------------------------------- */
265   count = 0;
266   for (K = firstPostorder(T); K != -1; K = nextPostorder(T, K))
267      {
268        vertex_root = first[K];
269        if (vertex_root == -1)
270          {
271            /* Should never happen */
272 #          if defined(PORD_INTSIZE64) || defined(INTSIZE64)
273            printf(" Internal error in mumps_pord, %ld\n",K);
274 #          else
275            printf(" Internal error in mumps_pord, %d\n",K);
276 #          endif
277            exit(-1);
278          }
279          /* for the principal column of the supervariable */
280        if (parent[K] == -1)
281          xadj_pe[vertex_root] = 0; /* root of the tree */
282        else
283          xadj_pe[vertex_root] = - (first[parent[K]]+1);
284        nv[vertex_root] = ncolfactor[K] + ncolupdate[K];
285        count++;
286        for (vertex = link[vertex_root]; vertex != -1; vertex = link[vertex])
287          /* for the secondary columns of the supervariable */
288          {
289            xadj_pe[vertex] = - (vertex_root+1);
290            nv[vertex] = 0;
291            count++;
292          }
293      }
294   /* ----------------------
295      free memory and return
296      ---------------------- */
297   free(first); free(link);
298   free(G->vwght);
299   free(G);
300   freeElimTree(T);
301   return (0);
302 }
303 #endif /* pord */
304