← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/examples/poisson_test/poisson_test.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/examples/poisson_test/poisson_test.c
 
1
/* zuse:    cc -o poisson_test -fast -native poisson_test.c mmio.c pcg.c -lm -xlic_lib=sunperf */
 
2
/* ru-lt13: cc -O poisson_test.c mmio.c pcg.c -o poisson_test -lblas -lg2c -lm */
 
3
#include <assert.h>
 
4
#include <stdio.h>
 
5
#include "mmio.h"
 
6
 
 
7
/* matrix A */
 
8
static int n_s;
 
9
static double *va_s, *da_s;
 
10
static int *ja_s, *ia_s;
 
11
 
 
12
/* CONVERT_COO_SSS - convert sparse matrix from COO to SSS format
 
13
 */
 
14
void convert_COO_SSS(int n, int nz,
 
15
                     int *i_coo, int *j_coo, double *v_coo,
 
16
                     int **ia, int **ja, double **va, double **da) {
 
17
  int i, k, l, t, nnz;
 
18
  int *root;
 
19
 
 
20
  root = (int *)malloc(n * sizeof(int));
 
21
  assert(root);
 
22
 
 
23
  /* allocate SSS matrix structure (1st part) */
 
24
  (*da) = (double *)malloc(n * sizeof(double));
 
25
 
 
26
  for (i = 0; i < n; i ++) {
 
27
    root[i] = 0;
 
28
    (*da)[i] = 0.0;
 
29
  }
 
30
 
 
31
  /* build n linked lists */
 
32
  nnz = 0;
 
33
  for (k = 0; k < nz; k ++) {
 
34
    if (i_coo[k] == j_coo[k])
 
35
      /* diagonal element */
 
36
      (*da)[i_coo[k]] = v_coo[k];
 
37
    else {
 
38
      /* off diagonal element */
 
39
      if (i_coo[k] < j_coo[k]) {   /* move to lower triangle */
 
40
        t = i_coo[k];
 
41
        i_coo[k] = j_coo[k];
 
42
        j_coo[k] = t;
 
43
      }
 
44
      i = i_coo[k];                /* link */
 
45
      i_coo[k] = root[i];
 
46
      root[i] = k;
 
47
      nnz ++;
 
48
    }
 
49
  }
 
50
 
 
51
  /* allocate SSS matrix structure (2nd part) */
 
52
  (*ia) = (int *)malloc((n+1) * sizeof(int));
 
53
  (*va) = (double *)malloc((nnz * sizeof(double)));
 
54
  (*ja) = (int *)malloc(nnz * sizeof(int));
 
55
 
 
56
  /* fill SSS matrix structure */
 
57
  k = 0;
 
58
  for (i = 0; i < n; i ++) {
 
59
    (*ia)[i] = k;
 
60
    l = root[i];
 
61
    while (l != 0) {
 
62
      (*ja)[k] = j_coo[l];
 
63
      (*va)[k] = v_coo[l];
 
64
      k ++;
 
65
      l = i_coo[l];
 
66
    }
 
67
  }
 
68
  (*ia)[n] = k;
 
69
  assert(k == nnz);
 
70
  free(root);
 
71
}
 
72
 
 
73
/* READ_MTX - read symmetric sparse matrix in MatrixMarket format
 
74
 */
 
75
void read_MTX_SSS(char *fname, int *n,
 
76
                  double **va, double **da, int **ja, int **ia) {
 
77
  int m, nz, ret_code, i;
 
78
  double *v_coo;
 
79
  int *i_coo, *j_coo;
 
80
  MM_typecode matcode;
 
81
  FILE *f;
 
82
 
 
83
  f = fopen(fname, "r");
 
84
  assert(f != NULL);
 
85
  ret_code = mm_read_banner(f, &matcode);
 
86
  assert(ret_code == 0);
 
87
  assert(mm_is_real(matcode) && mm_is_matrix(matcode) &&
 
88
         mm_is_sparse(matcode) && mm_is_symmetric(matcode));
 
89
  ret_code = mm_read_mtx_crd_size(f, &m, n, &nz);
 
90
  assert(ret_code == 0);
 
91
  assert(m == *n);
 
92
  /* read COO format */
 
93
  i_coo = (int *)malloc(nz * sizeof(int));
 
94
  j_coo = (int *)malloc(nz * sizeof(int));
 
95
  v_coo = (double *)malloc(nz * sizeof(double));
 
96
  assert(i_coo && j_coo && v_coo);
 
97
  for (i = 0; i < nz; i ++) {
 
98
    fscanf(f, "%d %d %lg\n", &i_coo[i], &j_coo[i], &v_coo[i]);
 
99
    i_coo[i]--;  /* adjust from 1-based to 0-based */
 
100
    j_coo[i]--;
 
101
  }
 
102
  fclose(f);
 
103
  /* convert to SSS format */
 
104
  convert_COO_SSS(*n, nz, i_coo, j_coo, v_coo, ia, ja, va, da);
 
105
  free(i_coo); free(j_coo); free(v_coo);
 
106
}
 
107
 
 
108
/* MATVEC - matrix vector multiplications
 
109
 */
 
110
void matvec(double *x, double *y) {
 
111
  double s, v, xi;
 
112
  int i, j, k;
 
113
 
 
114
  for (i = 0; i < n_s; i ++) {
 
115
    xi = x[i];
 
116
    s = 0.0;
 
117
    for (k = ia_s[i]; k < ia_s[i+1]; k ++) {
 
118
      j = ja_s[k];
 
119
      v = va_s[k];
 
120
      s += v * x[j];
 
121
      y[j] += v * xi;
 
122
    }
 
123
    y[i] = s + da_s[i]*xi;
 
124
  }
 
125
}
 
126
 
 
127
void main () {
 
128
  double *x, *b, *work;
 
129
  int i;
 
130
  double relres;
 
131
  int iter, flag;
 
132
 
 
133
  read_MTX_SSS("matrices/poi2d_100.mtx", &n_s, &va_s, &da_s, &ja_s, &ia_s);
 
134
 
 
135
  x = (double *) malloc(n_s * sizeof(double));
 
136
  b = (double *) malloc(n_s * sizeof(double));
 
137
  work = (double *) malloc(4*n_s * sizeof(double));
 
138
  assert(x != NULL && b != NULL && work != NULL);
 
139
 
 
140
  for (i = 0; i < n_s; i ++) {
 
141
    x[i] = 0.0;
 
142
    b[i] = 1.0;
 
143
  }
 
144
 
 
145
  printf("Starting PCG solver...\n");
 
146
  pcg(n_s, x, b, 1e-12, 2000, 1, &iter, &relres, &flag, work, matvec, NULL);
 
147
 
 
148
}