The R Project SVN R-packages

Rev

Rev 769 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 769 Rev 1381
1 
#include <Rdefines.h>
 1 
#include <Rdefines.h>
2 
#include "metis.h"
 2 
#include "metis.h"
3 
#include "Metis_utils.h"
 3 
#include "Metis_utils.h"
4 
 
 4 
 
- 
 
 5 
 
- 
 
 6 
 
- 
 
 7 
/**
- 
 
 8 
 * Establish a fill-reducing permutation for the sparse symmetric
- 
 
 9 
 * matrix of order n represented by the column pointers Tp and row
- 
 
 10 
 * indices Ti.
- 
 
 11 
 *
- 
 
 12 
 * @param n  order of the sparse symmetric matrix
- 
 
 13 
 * @param Tp  column pointers (total length n + 1)
- 
 
 14 
 * @param Ti  row indices (total length Tp[n])
- 
 
 15 
 * @param Perm array of length n to hold the permutation
- 
 
 16 
 * @param iPerm array of length n to hold the inverse permutation
- 
 
 17 
 *
- 
 
 18 
 */
5 
void ssc_metis_order(int n, const int Tp [], const int Ti [],
 19 
void ssc_metis_order(int n, const int Tp [], const int Ti [],
6 
		     int Perm[], int iPerm[])
 20 
		     int Perm[], int iPerm[])
7 
{
 21 
{
8 
    int  j, num_flag = 0, options_flag = 0;
 22 
    int  j, num_flag = 0, options_flag = 0;
9 
    idxtype
 23 
    idxtype
10 
	*perm = Calloc(n, idxtype), /* in case idxtype != int */
 24 
	*perm = Calloc(n, idxtype), /* in case idxtype != int */
11 
	*iperm = Calloc(n, idxtype),
 25 
	*iperm = Calloc(n, idxtype),
12 
	*xadj = Calloc(n+1, idxtype),
 26 
	*xadj = Calloc(n+1, idxtype),
13 
	*adj = Calloc(2 * (Tp[n] - n), idxtype);
 27 
	*adj = Calloc(2 * (Tp[n] - n), idxtype);
14 
 
 28 
 
- 
 
 29 
				/* check row indices for correct range */
- 
 
 30 
    for (j = 0; j < Tp[n]; j++)
- 
 
 31 
      if (Ti[j] < 0 || Ti[j] >= n)
- 
 
 32 
	error(_("row index Ti[%d] = %d is out of range [0,%d]"),
- 
 
 33 
	      j, Ti[j], n - 1);
15 
				/* temporarily use perm to store lengths */
 34 
				/* temporarily use perm to store lengths */
16 
    memset(perm, 0, sizeof(idxtype) * n);
 35 
    AZERO(perm, n);
17 
    for (j = 0; j < n; j++) {
 36 
    for (j = 0; j < n; j++) {
18 
	int ip, p2 = Tp[j+1];
 37 
	int ip, p2 = Tp[j+1];
19 
	for (ip = Tp[j]; ip < p2; ip++) {
 38 
	for (ip = Tp[j]; ip < p2; ip++) {
20 
	    int i = Ti[ip];
 39 
	    int i = Ti[ip];
21 
	    if (i != j) {
 40 
	    if (i != j) {
22 
		perm[i]++;
 41 
		perm[i]++;
23 
		perm[j]++;
 42 
		perm[j]++;
24 
	    }
 43 
	    }
25 
	}
 44 
	}
26 
    }
 45 
    }
27 
    xadj[0] = 0;
 46 
    xadj[0] = 0;
28 
    for (j = 0; j < n; j++) xadj[j+1] = xadj[j] + perm[j];
 47 
    for (j = 0; j < n; j++) xadj[j+1] = xadj[j] + perm[j];
29 
				/* temporarily use perm to store pointers */
 48 
				/* temporarily use perm to store pointers */
30 
    Memcpy(perm, xadj, n);
 49 
    Memcpy(perm, xadj, n);
31 
    for (j = 0; j < n; j++) {
 50 
    for (j = 0; j < n; j++) {
32 
	int ip, p2 = Tp[j+1];
 51 
	int ip, p2 = Tp[j+1];
33 
	for (ip = Tp[j]; ip < p2; ip++) {
 52 
	for (ip = Tp[j]; ip < p2; ip++) {
34 
	    int i = Ti[ip];
 53 
	    int i = Ti[ip];
35 
	    if (i != j) {
 54 
	    if (i != j) {
36 
		adj[perm[i]] = j;
 55 
		adj[perm[i]] = j;
37 
		adj[perm[j]] = i;
 56 
		adj[perm[j]] = i;
38 
		perm[i]++;
 57 
		perm[i]++;
39 
		perm[j]++;
 58 
		perm[j]++;
40 
	    }
 59 
	    }
41 
	}
 60 
	}
42 
    }
 61 
    }
43 
    METIS_NodeND(&n, xadj, adj, &num_flag, &options_flag, perm, iperm);
 62 
    METIS_NodeND(&n, xadj, adj, &num_flag, &options_flag, perm, iperm);
44 
    for (j = 0; j < n; j++) {
 63 
    for (j = 0; j < n; j++) {
45 
	Perm[j] = (int) perm[j];
 64 
	Perm[j] = (int) perm[j];
46 
	iPerm[j] = (int) iperm[j];
 65 
	iPerm[j] = (int) iperm[j];
47 
    }
 66 
    }
48 
    Free(iperm); Free(perm); Free(xadj); Free(adj);
 67 
    Free(iperm); Free(perm); Free(xadj); Free(adj);
49 
}
 68 
}