*************************************************************************/
#define STANDALONE
-#include "SparseMatrix.h"
#include "country_graph_coloring.h"
#include "power.h"
+/* #include "general.h" */
+/* #include "SparseMatrix.h" */
+/* #include "country_graph_coloring.h" */
+/* #include "power.h" */
+#include "PriorityQueue.h"
+#include "time.h"
+
+/* int Verbose = FALSE; */
+
static void get_local_12_norm(int n, int i, int *ia, int *ja, int *p, real *norm){
- int j;
+ int j, nz;
norm[0] = n; norm[1] = 0;
for (j = ia[i]; j < ia[i+1]; j++){
if (ja[j] == i) continue;
norm[0] = MIN(norm[0], ABS(p[i] - p[ja[j]]));
+ nz++;
norm[1] += ABS(p[i] - p[ja[j]]);
}
+ if (nz > 0) norm[1] /= nz;
}
static void get_12_norm(int n, int *ia, int *ja, int *p, real *norm){
- int i, j;
- norm[0] = n; norm[1] = 0;
+ /* norm[0] := antibandwidth
+ norm[1] := (\sum_{{i,j}\in E} |p[i] - p[j]|)/|E|
+ norm[2] := (\sum_{i\in V} (Min_{{j,i}\in E} |p[i] - p[j]|)/|V|
+ */
+ int i, j, nz = 0;
+ real tmp;
+ norm[0] = n; norm[1] = 0; norm[2] = 0;
for (i = 0; i < n; i++){
+ tmp = n;
for (j = ia[i]; j < ia[i+1]; j++){
if (ja[j] == i) continue;
norm[0] = MIN(norm[0], ABS(p[i] - p[ja[j]]));
norm[1] += ABS(p[i] - p[ja[j]]);
+ tmp = MIN(tmp, ABS(p[i] - p[ja[j]]));
+ nz++;
}
+ norm[2] += tmp;
}
+ norm[2] /= n;
+ norm[1] /= nz;
}
-void country_graph_coloring(int seed, SparseMatrix A, int **p, real *norm_1){
+static void update_pmin_pmax_aband(int n, int u, int *ia, int *ja, int *p, int *pmin, int *pmax, int *aband_local){
+ int j, aband_u;
+ pmin[u] = n; pmax[u] = -1; aband_u = n;
+ for (j = ia[u]; j < ia[u+1]; j++) {
+ if (ja[j] == u) continue;
+ pmin[u] = MIN(pmin[u], ABS(p[u] - p[ja[j]]));
+ pmax[u] = MIN(pmax[u], ABS(p[u] - p[ja[j]]));
+ aband_u = MIN(aband_u, ABS(p[u] - p[ja[j]]));
+ }
+ aband_local[u] = aband_u;
+}
+
+
+static int check_swap(int n, int *ia, int *ja,
+ int u, int p_u, int v, int p_v, int *aband_local, int *p, int *p_inv, int aband, PriorityQueue pq, int *pmax, int *pmin, real lambda){
+ /* check if u should swap with v to improve u, without demaging v. Return TRUE if swap is successful. FALSE otherwise. */
+
+ int j, aband_v, aband_v1, aband_u, aband_u1;
+
+ aband_v = aband_local[v];
+ aband_u = aband_local[u];
+
+ /* if swaping u and v makes v worse & becomes/remains critical, don't do. We first quick check using the max/min neighbor indices.
+ No need to check the other way around since the calling function have ensured that.
+ */
+ if (ABS(p_u - pmin[v]) < aband_v && ABS(p_u - pmin[v]) <= lambda*aband) return FALSE;
+ if (ABS(p_u - pmax[v]) < aband_v && ABS(p_u - pmax[v]) <= lambda*aband) return FALSE;
+
+ /* now check in details whether v should swap to u. Do not accept if this makes the antiband width of u worse */
+ aband_u1 = n;
+ for (j = ja[u]; j < ja[u+1]; j++){
+ if (ja[j] == u) continue;
+ if (ABS(p_v - p[ja[j]]) < aband_u) {
+ return FALSE;
+ }
+ aband_u1 = MIN(aband_u1, ABS(p_v - p[ja[j]]));
+ }
+
+ /* now check in details whether u should swap to v. Do not accept if this makes antibandwidth of v worse && make/keep v in the critical group */
+ aband_v1 = n;
+ for (j = ja[v]; j < ja[v+1]; j++){
+ if (ja[j] == v) continue;
+ if (ABS(p_u - p[ja[j]]) < aband_v && ABS(p_u - p[ja[j]]) <= lambda*aband) {
+ return FALSE;
+ }
+ aband_v1 = MIN(aband_v1, ABS(p_u - p[ja[j]]));
+ }
+
+ /* now check if local antiband width has been improved. By that we mean u is improved, or u unchanged, but v improved. */
+ assert(aband_u1 >= aband_u);
+ if (aband_u1 > aband_u || (aband_u1 == aband_u && aband_v1 > aband_v)){
+ p[u] = p_v;
+ p[v] = p[u];
+ p_inv[p[u]] = u;
+ p_inv[p[v]] = v;
+
+ update_pmin_pmax_aband(n, u, ia, ja, p, pmin, pmax, aband_local);
+ update_pmin_pmax_aband(n, v, ia, ja, p, pmin, pmax, aband_local);
+
+ /* this may be expensive, but I see no way to keep pmin/pmax/aband_local consistent other than this update of u/v's neighbors */
+ for (j = ia[u]; j < ia[u+1]; j++) {
+ update_pmin_pmax_aband(n, ja[j], ia, ja, p, pmin, pmax, aband_local);
+ }
+
+ for (j = ia[u]; j < ia[u+1]; j++) {
+ update_pmin_pmax_aband(n, ja[j], ia, ja, p, pmin, pmax, aband_local);
+ }
+ return TRUE;
+ }
+
+
+ return FALSE;
+}
+
+void improve_antibandwidth_by_swapping_cheap(SparseMatrix A, int *p){
+ /*on entry:
+ A: the graph, must be symmetric matrix
+ p: a permutation array of length n
+ lambda: threshold for deciding critical vertices.*/
+ real lambda = 1.2;
+ int n = A->m;
+ int *p_inv;
+ int i, j;
+ int *pmax, *pmin;/* max and min index of neighbors */
+ int *ia = A->ia, *ja = A->ja;
+ int aband = n;/* global antibandwidth*/
+ int *aband_local;/* antibandwidth for each node */
+ PriorityQueue pq = NULL;
+ int progress = TRUE, u, v, gain, aband_u, p_u, p_v, swapped;
+
+ pq = PriorityQueue_new(n, n);
+ p_inv = MALLOC(sizeof(int)*n);
+ pmax = MALLOC(sizeof(int)*n);
+ pmin = MALLOC(sizeof(int)*n);
+ aband_local = MALLOC(sizeof(int)*n);
+
+ while (progress) {
+ progress = FALSE;
+ for (i = 0; i < n; i++){
+ pmax[i] = -1; pmin[i] = n+1;
+ assert(p[i] >= 0 && p[i] < n);
+ p_inv[p[i]] = i;
+ aband_local[i] = n;
+ for (j = ia[i]; j < ia[i+1]; j++){
+ if (ja[j] == i) continue;
+ pmax[i] = MAX(pmax[i], p[ja[j]]);
+ pmin[i] = MIN(pmin[i], p[ja[j]]);
+ aband_local[i] = MIN(aband_local[i], ABS(p[i] - p[ja[j]]));
+ }
+ aband = MIN(aband, aband_local[i]);
+ }
+ fprintf(stderr," antibandwidth = %d", aband);
+
+ /* find critical vertices */
+ for (i = 0; i < n; i++){
+ if (aband_local[i] <= lambda*aband){
+ PriorityQueue_push(pq, i, n - aband_local[i]);
+ }
+ }
+
+ /* check critcal nodes u to see if any swaps with v are possible */
+ while (PriorityQueue_pop(pq, &u, &gain)){
+ aband_u = n - gain;
+ p_u = p[u];
+ assert(aband_u <= lambda*aband);
+ assert(aband_u == aband_local[u]);
+ swapped = FALSE;
+ for (p_v = 0; p_v <= pmin[u] - aband_u; p_v++){
+ v = p_inv[p_v];
+ if (check_swap(n, ia, ja, u, p_u, v, p_v, aband_local, p, p_inv, aband, pq, pmax, pmin, lambda)){
+ swapped = TRUE; progress = TRUE;
+ break;
+ }
+ }
+ if (swapped) continue;
+
+ for (p_v = pmax[u] + aband_u; p_v < n; p_v++){
+ v = p_inv[p_v];
+ if (check_swap(n, ia, ja, u, p_u, v, p_v, aband_local, p, p_inv, aband, pq, pmax, pmin, lambda)){
+ swapped = TRUE; progress = TRUE;
+ break;
+ }
+ }
+ if (swapped) continue;
+
+ for (p_v = pmin[u] + aband_u; p_v <= pmax[u] - aband_u; p_v++) {
+ v = p_inv[p_v];
+ if (check_swap(n, ia, ja, u, p_u, v, p_v, aband_local, p, p_inv, aband, pq, pmax, pmin, lambda)){
+ progress = TRUE;
+ break;
+ }
+ }
+
+
+ }
+ }
+
+
+ FREE(p_inv);
+ FREE(pmax);
+ FREE(pmin);
+ FREE(aband_local);
+ PriorityQueue_delete(pq);
+
+}
+
+void improve_antibandwidth_by_swapping(SparseMatrix A, int *p){
+ int improved = TRUE;
+ int n = A->m, i, j, *ia = A->ia, *ja = A->ja;
+ real norm = n, norm1[3], norm2[3], norm11[3], norm22[3];
+ real pi, pj;
+ real start = clock();
+ FILE *fp;
+
+ if (Verbose){
+ fprintf(stderr,"saving timing vs antiband data to timing_greedy\n");
+ fp = fopen("timing_greedy","w");
+ }
+ assert(SparseMatrix_is_symmetric(A, TRUE));
+ while (improved){
+ improved = FALSE; norm = n;
+ for (i = 0; i < n; i++){
+ get_local_12_norm(n, i, ia, ja, p, norm1);
+ for (j = 0; j < n; j++){
+ if (j == i) continue;
+ get_local_12_norm(n, j, ia, ja, p, norm2);
+ norm = MIN(norm, norm2[0]);
+ pi = (p)[i]; pj = (p)[j];
+ (p)[i] = pj;
+ (p)[j] = pi;
+ get_local_12_norm(n, i, ia, ja, p, norm11);
+ get_local_12_norm(n, j, ia, ja, p, norm22);
+ if (MIN(norm11[0],norm22[0]) > MIN(norm1[0],norm2[0])){
+ // ||
+ //(MIN(norm11[0],norm22[0]) == MIN(norm1[0],norm2[0]) && norm11[1]+norm22[1] > norm1[1]+norm2[1])) {
+ improved = TRUE;
+ norm1[0] = norm11[0];
+ norm1[1] = norm11[1];
+ continue;
+ }
+ (p)[i] = pi;
+ (p)[j] = pj;
+ }
+ if (i%100 == 0 && Verbose) {
+ get_12_norm(n, ia, ja, p, norm1);
+ fprintf(fp,"%f %f %f\n", (real) (clock() - start)/(CLOCKS_PER_SEC), norm1[0], norm1[2]);
+ }
+ }
+ if (Verbose) {
+ get_12_norm(n, ia, ja, p, norm1);
+ fprintf(stderr, "aband = %f, aband_avg = %f\n", norm1[0], norm1[2]);
+ fprintf(fp,"%f %f %f\n", (real) (clock() - start)/(CLOCKS_PER_SEC), norm1[0], norm1[2]);
+ }
+ }
+}
+
+void country_graph_coloring_internal(int seed, SparseMatrix A, int **p, real *norm_1, int do_swapping){
int n = A->m, i, j, jj;
SparseMatrix L, A2;
int *ia = A->ia, *ja = A->ja;
int a = -1.;
real nrow;
real *v = NULL;
- real norm1[2], norm2[2], norm11[2], norm22[2], norm = n;
- real pi, pj;
- int improved = TRUE;
-
+ real norm1[3];
+ clock_t start, start2;
+ start = clock();
assert(A->m == A->n);
A2 = SparseMatrix_symmetrize(A, TRUE);
ia = A2->ia; ja = A2->ja;
{
int maxit = 100;
real tol = 0.00001;
- power_method(L, seed, maxit, tol, &v);
+ real eig, *eigv;
+ eigv = &eig;
+ power_method(matvec_sparse, L, L->n, 1, seed, maxit, tol, &v, &eigv);
}
vector_ordering(n, v, p, TRUE);
+ fprintf(stderr, "cpu time for spectral ordering (before greedy) = %f\n", (real) (clock() - start)/(CLOCKS_PER_SEC));
+ start2 = clock();
/* swapping */
- while (improved){
- improved = FALSE; norm = n;
- for (i = 0; i < n; i++){
- get_local_12_norm(n, i, ia, ja, *p, norm1);
- for (j = 0; j < n; j++){
- if (j == i) continue;
- get_local_12_norm(n, j, ia, ja, *p, norm2);
- norm = MIN(norm, norm2[0]);
- pi = (*p)[i]; pj = (*p)[j];
- (*p)[i] = pj;
- (*p)[j] = pi;
- get_local_12_norm(n, i, ia, ja, *p, norm11);
- get_local_12_norm(n, j, ia, ja, *p, norm22);
- if (MIN(norm11[0],norm22[0]) > MIN(norm1[0],norm2[0])){
- // ||
- //(MIN(norm11[0],norm22[0]) == MIN(norm1[0],norm2[0]) && norm11[1]+norm22[1] > norm1[1]+norm2[1])) {
- improved = TRUE;
- norm1[0] = norm11[0];
- norm1[1] = norm11[1];
- continue;
- }
- (*p)[i] = pi;
- (*p)[j] = pj;
- }
- }
- if (Verbose) {
- get_12_norm(n, ia, ja, *p, norm1);
- fprintf(stderr, "norm = %f", norm);
- fprintf(stderr, "norm1 = %f, norm2 = %f\n", norm1[0], norm1[1]);
+ if (do_swapping) {
+ if (do_swapping == DO_SWAPPING){
+ improve_antibandwidth_by_swapping(A2, *p);
+ } else if (do_swapping == DO_SWAPPING_CHEAP) {
+ improve_antibandwidth_by_swapping_cheap(A2, *p);
+ } else {
+ assert(0);
}
}
+ fprintf(stderr, "cpu time for greedy refinement = %f\n", (real) (clock() - start2)/(CLOCKS_PER_SEC));
+
+ fprintf(stderr, "cpu time for spectral + greedy = %f\n", (real) (clock() - start)/(CLOCKS_PER_SEC));
+
get_12_norm(n, ia, ja, *p, norm1);
*norm_1 = norm1[0];
if (A2 != A) SparseMatrix_delete(A2);
SparseMatrix_delete(L);
}
+void country_graph_coloring(int seed, SparseMatrix A, int **p, real *norm_1){
+ country_graph_coloring_internal(seed, A, p, norm_1, DO_SWAPPING);
+}
+
+void country_graph_coloring_general(int seed, SparseMatrix A, int **p, real *norm_1, int do_swapping){
+ country_graph_coloring_internal(seed, A, p, norm_1, do_swapping);
+}
+
+
+void improve_antibandwidth_by_swapping_for_fortran(int *n, int *nz, int *ja, int *ia, int *p, int *aprof, int *verbose){
+ /* n: dimension
+ nz: number of nonzeros
+ ja: array of size nz, holds column indices. 1- based
+ ia: array of size n+1, position of row starts. 1-based.
+ p: holds the position of row/col i in the permuted matrix. 1-based
+ aprof: aprof[0] is the antiband width on entry. aprof[1] is the abtiband on exit
+ Verbose: a flag, when set to nonzero, will print the input matrix, as well as others things.
+ */
+ SparseMatrix A, A2;
+ real norm1[3];
+ int i, j, jj;
+ clock_t start;
+ real cpu;
+
+ Verbose = *verbose;
+ A = SparseMatrix_new(*n, *n, 1, MATRIX_TYPE_PATTERN, FORMAT_COORD);
+
+ for (i = 0; i < *n; i++){
+ for (j = ia[i] - 1; j < ia[i+1] - 1; j++){
+ jj = ja[j] - 1;
+ A = SparseMatrix_coordinate_form_add_entries(A, 1, &i, &jj, NULL);
+ }
+ }
+ A2 = SparseMatrix_from_coordinate_format(A);
+ if (Verbose && 0) SparseMatrix_print("A = ",A2);
+
+ SparseMatrix_delete(A);
+ A = SparseMatrix_symmetrize(A2, TRUE);
+
+ for (i = 0; i < *n; i++) p[i]--;
+
+ ia = A->ia; ja = A->ja;
+ get_12_norm(*n, ia, ja, p, norm1);
+ if (Verbose) fprintf(stderr,"on entry antibandwidth = %f\n", norm1[0]);
+ aprof[0] = (int) norm1[0];
+
+ start = clock();
+ improve_antibandwidth_by_swapping(A, p);
+ cpu = (clock() - start)/(CLOCKS_PER_SEC);
+ fprintf(stderr,"cpu = %f\n", cpu);
+
+ get_12_norm(*n, ia, ja, p, norm1);
+ if (Verbose) fprintf(stderr,"on exit antibandwidth = %f\n", norm1[0]);
+ aprof[1] = (int) norm1[0];
+ SparseMatrix_delete(A);
+ SparseMatrix_delete(A2);
+ for (i = 0; i < *n; i++) p[i]++;
+}
+
+void improve_antibandwidth_by_swapping_for_fortran_(int *n, int *nz, int *ja, int *ia, int *p, int *aprof, int *Verbose){
+ improve_antibandwidth_by_swapping_for_fortran(n, nz, ja, ia, p, aprof, Verbose);
+}
+void IMPROVE_ANTIBANDWIDTH_BY_SWAPPING_FOR_FORTRAN(int *n, int *nz, int *ja, int *ia, int *p, int *aprof, int *Verbose){
+ improve_antibandwidth_by_swapping_for_fortran(n, nz, ja, ia, p, aprof, Verbose);
+}
+void IMPROVE_ANTIBANDWIDTH_BY_SWAPPING_FOR_FORTRAN_(int *n, int *nz, int *ja, int *ia, int *p, int *aprof, int *Verbose){
+ improve_antibandwidth_by_swapping_for_fortran(n, nz, ja, ia, p, aprof, Verbose);
+}
+
* Contributors: See CVS logs. Details at http://www.graphviz.org/
*************************************************************************/
+#include "power.h"
#include "SparseMatrix.h"
-/* #include "matrix_market.h" */
-void power_method(SparseMatrix A, int random_seed, int maxit, real tol, real **eigv){
- /* find the largest eigenvector of a matrix A. Result in eigv. if eigv == NULL; memory will be allocated.
+
+void power_method(void (*matvec)(void *, int, int, real*, real **, int, int*),
+ void *A, int n, int K, int random_seed, int maxit, real tol, real **eigv, real **eigs){
+ /* find k-largest eigenvectors of a matrix A. Result in eigv. if eigv == NULL; memory will be allocated.
maxium of maxit iterations will be done, and tol is the convergence criterion
- This converges only if the largest eigenvector/value is real and the next largest eigenvalue separate from teh largest one
+ This converges only if the largest eigenvectors/values are real (e.g., if A is symmetric) and the
+ next largest eigenvalues separate from the largest ones
+
+ input:
+ matvec: a function point that takes a matrix M and a vector u, produce v = M.u
+ A: the matrix
+ n: dimension of matrix A
+ K: number of eigenes to find
+ random_seed: seed for eigenvector initialization
+ matrix: max number f iterations
+ tol: accuracy control
+
+ output:
+ eigv: eigenvectors. The i-th is at eigvs[i*n, i*(n+1) - 1]
+ eigs: eigenvalues. The i-th is at eigs[i]
+
+ Function PowerIteration (A – m × m matrix )
+ % This function computes u1, u2, . . . , uk, the first k eigenvectors of S.
+ const tol ← 0.001
+ for i = 1 to k do
+ . ui ← random
+ . ui ← ui/||ui||
+ . do
+ . vi ← ui
+ . % orthogonalize against previous eigenvectors
+ . for j = 1 to i − 1 do
+ . vi ← vi − (vi^Tvi)vj
+ . end for
+ . ui ← A vi/||A vi||
+ . while (ui^T vi < 1-tol) (halt when direction change is small)
+ . vi = ui
+ end for
+ return v1,v2,...
*/
- int n;
- real *v, *u;
+ real **v, *u, *vv;
int iter = 0;
real res, unorm;
- int i;
- assert(A->m == A->n);
- assert(A->type = MATRIX_TYPE_REAL || A->type == MATRIX_TYPE_INTEGER);
+ int i, j, k;
+ real uij;
+ int flag;
- n = A->m;
- if (!(*eigv)) *eigv = MALLOC(sizeof(real)*n);
- u = MALLOC(sizeof(real)*n);
+ K = MAX(0, MIN(n, K));
+ assert(K <= n && K > 0);
- srand(random_seed);
+ if (!(*eigv)) *eigv = MALLOC(sizeof(real)*n*K);
+ if (!(*eigs)) *eigs = MALLOC(sizeof(real)*K);
+ v = MALLOC(sizeof(real*)*K);
- for (i = 0; i < n; i++) (*eigv)[i] = drand();
+ vv = MALLOC(sizeof(real)*n);
+ u = MALLOC(sizeof(real)*n);
- res = vector_product(n, *eigv, *eigv);
- if (res > 0) res = 1/res;
- for (i = 0; i < n; i++) (*eigv)[i] = (*eigv)[i]*res;
-
- v = *eigv;
-
- do {
- SparseMatrix_multiply_vector(A, v, &u, FALSE);
+ srand(random_seed);
- unorm = vector_product(n, u, u);/* ||u||^2 */
- unorm = sqrt(unorm);
- if (unorm > 0) unorm = 1/unorm;
- res = 0.;
+ for (k = 0; k < K; k++){
+ //fprintf(stderr,"calculating eig k ==================== %d\n",k);
+ v[k] = &((*eigv)[k*n]);
+ for (i = 0; i < n; i++) u[i] = drand();
+ res = sqrt(vector_product(n, u, u));
+ if (res > 0) res = 1/res;
for (i = 0; i < n; i++) {
- u[i] = u[i]*unorm;
- res = res + (u[i] - v[i])*(u[i] - v[i]);
- v[i] = u[i];
+ u[i] = u[i]*res;
+ v[k][i] = u[i];
}
-
- /*
- printf("=== %d === %f\n",iter,res);
- printf("{");
- {int j;
- for (j = 0; j < MIN(10,n); j++){
- if (j == n-1){
- printf("%f",v[j]);
+ /*
+ fprintf(stderr,"inital vec=");
+ for (i = 0; i < n; i++) fprintf(stderr,"%f,",u[i]);fprintf(stderr,"\n");
+ */
+ iter = 0;
+ do {
+
+
+ /* normalize against previous eigens */
+ for (j = 0; j < k; j++){
+ uij = vector_product(n, u, v[j]);
+ for (i = 0; i < n; i++) {
+ u[i] = u[i] - uij *v[j][i];
+ }
+ }
+ matvec(A, n, n, u, &vv, FALSE, &flag);
+ assert(!flag);
+
+ /*
+ fprintf(stderr,"normalized aginst prev vec=");
+ for (i = 0; i < n; i++) fprintf(stderr,"%f,",u[i]);fprintf(stderr,"\n");
+ */
+
+ unorm = vector_product(n, vv, vv);/* ||u||^2 */
+ unorm = sqrt(unorm);
+ (*eigs)[k] = unorm;
+ if (unorm > 0) {
+ unorm = 1/unorm;
} else {
- printf("%f,",v[j]);
+ // ||A.v||=0, so v must be an eigenvec correspond to eigenvalue zero
+ for (i = 0; i < n; i++) vv[i] = u[i];
+ unorm = sqrt(vector_product(n, vv, vv));
+ if (unorm > 0) unorm = 1/unorm;
}
- }
- printf("\n");
- }
- */
+ res = 0.;
- } while (res/n > tol*tol && iter++ < maxit);
-
-
+ for (i = 0; i < n; i++) {
+ //res = MAX(res, ABS(vv[i]-(*eigs)[k]*u[i]));
+ u[i] = vv[i]*unorm;
+ res = res + u[i]*v[k][i];
+ v[k][i] = u[i];
+ }
+ //fprintf(stderr,"res=%g, tol = %g, res < 1-tol=%d\n",res, tol,res < 1 - tol);
+ } while (res < 1 - tol && iter++ < maxit);
+ //} while (iter++ < maxit);
+ //fprintf(stderr,"iter= %d, res=%f\n",iter, res);
+ }
+ FREE(u);
+ FREE(vv);
}
-/*
-main(){
- real *v = NULL;
- int i;
- int n;
+void matvec_sparse(void *M, int m, int n, real *u, real **v, int transpose,
+ int *flag){
SparseMatrix A;
- A = SparseMatrix_import_matrix_market(stdin, FORMAT_CSR);
- n = A->m;
-
- power_method(A, 123, 100, 0.00001, &v);
+ *flag = 0;
+ A = (SparseMatrix) M;
+ SparseMatrix_multiply_vector(A, u, v, transpose);
+ return;
+}
- printf("{");
- for (i = 0; i < n; i++){
- if (i == n-1){
- printf("%f",v[i]);
+void mat_print_dense(real *M, int m, int n){
+ int i, j;
+ fprintf(stderr,"{");
+ for (i = 0; i < m; i++){
+ fprintf(stderr,"{");
+ for (j = 0; j < n; j++){
+ if (j != 0) fprintf(stderr,",");
+ fprintf(stderr,"%f",M[i*n+j]);
+ }
+ if (i != m-1){
+ fprintf(stderr,"},\n");
} else {
- printf("%f,",v[i]);
+ fprintf(stderr,"}");
}
}
- printf("}\n");
+ fprintf(stderr,"}\n");
+}
+
+void matvec_dense(void *M, int m, int n, real *u, real **v, int transpose,
+ int *flag){
+ /* M.u or M^T.u */
+ real *A;
+ int i, j;
+
+ A = (real*) M;
+ *flag = 0;
+
+
+ if (!transpose){
+ if (!(*v)) *v = MALLOC(sizeof(real)*m);
+ for (i = 0; i < m; i++){
+ (*v)[i] = 0;
+ for (j = 0; j < n; j++){
+ (*v)[i] += A[i*n+j]*u[j];
+ }
+ }
+ } else {
+ if (!(*v)) *v = MALLOC(sizeof(real)*n);
+ for (i = 0; i < n; i++){
+ (*v)[i] = 0;
+ }
+ for (i = 0; i < m; i++){
+ for (j = 0; j < n; j++){
+ (*v)[j] += A[i*n+j]*u[i];
+ }
+ }
+ }
+
+ return;
}
-*/
projects / graphviz / commitdiff