The R Project SVN R-packages

Rev

Rev 769 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 769 Rev 1381
Line 1... Line 1...
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) {