matrix_ops.c

Go to the documentation of this file.

/*************************************************************************
 * Copyright (c) 2011 AT&T Intellectual Property
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v2.0
 * which accompanies this distribution, and is available at
 * https://www.eclipse.org/org/documents/epl-2.0/EPL-2.0.html
 *
 * Contributors: Details at https://graphviz.org
 *************************************************************************/
 
#include "config.h"
 
#include <math.h>
#include <neatogen/matrix_ops.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <util/alloc.h>
 
static const double p_iteration_threshold = 1e-3;
 
bool power_iteration(double *const *square_mat, int n, int neigs,
                     double **eigs) {
  /* compute the 'neigs' top eigenvectors of 'square_mat' using power iteration
   */
 
  int i;
  double *tmp_vec = gv_calloc(n, sizeof(double));
  double *last_vec = gv_calloc(n, sizeof(double));
  double angle;
  int iteration = 0;
  const int Max_iterations = 30 * n;
 
  double tol = 1 - p_iteration_threshold;
 
  if (neigs >= n) {
    neigs = n;
  }
 
  double *const evals = gv_calloc(neigs, sizeof(evals[0]));
  for (i = 0; i < neigs; i++) {
    double *const curr_vector = eigs[i];
    /* guess the i-th eigen vector */
  choose:
    for (int j = 0; j < n; j++)
      curr_vector[j] = rand() % 100;
    /* orthogonalize against higher eigenvectors */
    for (int j = 0; j < i; j++) {
      const double alpha = -vectors_inner_product(n, eigs[j], curr_vector);
      scadd(curr_vector, n - 1, alpha, eigs[j]);
    }
    double len = norm(curr_vector, n - 1);
    if (len < 1e-10) {
      // we have chosen a vector colinear with previous ones
      goto choose;
    }
    vectors_scalar_mult(n, curr_vector, 1.0 / len, curr_vector);
    iteration = 0;
    do {
      iteration++;
      copy_vector(n, curr_vector, last_vec);
 
      right_mult_with_vector_d(square_mat, n, n, curr_vector, tmp_vec);
      copy_vector(n, tmp_vec, curr_vector);
 
      /* orthogonalize against higher eigenvectors */
      for (int j = 0; j < i; j++) {
        const double alpha = -vectors_inner_product(n, eigs[j], curr_vector);
        scadd(curr_vector, n - 1, alpha, eigs[j]);
      }
      len = norm(curr_vector, n - 1);
      if (len < 1e-10 || iteration > Max_iterations) {
        /* We have reached the null space (e.vec. associated with e.val. 0) */
        goto exit;
      }
 
      vectors_scalar_mult(n, curr_vector, 1.0 / len, curr_vector);
      angle = vectors_inner_product(n, curr_vector, last_vec);
    } while (fabs(angle) < tol);
    evals[i] = angle * len; /* this is the Rayleigh quotient (up to errors due
                               to orthogonalization): u*(A*u)/||A*u||)*||A*u||,
                               where u=last_vec, and ||u||=1
                             */
  }
exit:
  for (; i < neigs; i++) {
    /* compute the smallest eigenvector, which are  */
    /* probably associated with eigenvalue 0 and for */
    /* which power-iteration is dangerous */
    double *const curr_vector = eigs[i];
    /* guess the i-th eigen vector */
    for (int j = 0; j < n; j++)
      curr_vector[j] = rand() % 100;
    /* orthogonalize against higher eigenvectors */
    for (int j = 0; j < i; j++) {
      const double alpha = -vectors_inner_product(n, eigs[j], curr_vector);
      scadd(curr_vector, n - 1, alpha, eigs[j]);
    }
    const double len = norm(curr_vector, n - 1);
    vectors_scalar_mult(n, curr_vector, 1.0 / len, curr_vector);
    evals[i] = 0;
  }
 
  /* sort vectors by their evals, for overcoming possible mis-convergence: */
  for (i = 0; i < neigs - 1; i++) {
    int largest_index = i;
    double largest_eval = evals[largest_index];
    for (int j = i + 1; j < neigs; j++) {
      if (largest_eval < evals[j]) {
        largest_index = j;
        largest_eval = evals[largest_index];
      }
    }
    if (largest_index != i) { /* exchange eigenvectors: */
      copy_vector(n, eigs[i], tmp_vec);
      copy_vector(n, eigs[largest_index], eigs[i]);
      copy_vector(n, tmp_vec, eigs[largest_index]);
 
      evals[largest_index] = evals[i];
      evals[i] = largest_eval;
    }
  }
 
  free(evals);
  free(tmp_vec);
  free(last_vec);
 
  return iteration <= Max_iterations;
}
 
void mult_dense_mat_d(double **A, float **B, int dim1, int dim2, int dim3,
                      double ***CC) {
  // A is dim1 × dim2, B is dim2 × dim3, C = A × B
 
  int i, j, k;
  double sum;
 
  double *storage = gv_calloc(dim1 * dim3, sizeof(double));
  double **C = *CC = gv_calloc(dim1, sizeof(double *));
 
  for (i = 0; i < dim1; i++) {
    C[i] = storage;
    storage += dim3;
  }
 
  for (i = 0; i < dim1; i++) {
    for (j = 0; j < dim3; j++) {
      sum = 0;
      for (k = 0; k < dim2; k++) {
        sum += A[i][k] * B[k][j];
      }
      C[i][j] = sum;
    }
  }
}
 
void mult_sparse_dense_mat_transpose(vtx_data *A, double **B, int dim1,
                                     int dim2, float ***CC) {
  // A is dim1 × dim1 and sparse, B is dim2 × dim1, C = A × B
 
  int i, j;
  double sum;
  float *ewgts;
  int *edges;
  float *storage = gv_calloc(dim1 * dim2, sizeof(A[0]));
  float **C = *CC = gv_calloc(dim1, sizeof(A));
 
  for (i = 0; i < dim1; i++) {
    C[i] = storage;
    storage += dim2;
  }
 
  for (i = 0; i < dim1; i++) {
    edges = A[i].edges;
    ewgts = A[i].ewgts;
    const size_t nedges = A[i].nedges;
    for (j = 0; j < dim2; j++) {
      sum = 0;
      for (size_t k = 0; k < nedges; k++) {
        sum += ewgts[k] * B[j][edges[k]];
      }
      C[i][j] = (float)(sum);
    }
  }
}
 
/* Scaled add - fills double vec1 with vec1 + alpha*vec2 over range*/
void scadd(double *vec1, int end, double fac, double *vec2) {
  int i;
 
  for (i = end + 1; i; i--) {
    (*vec1++) += fac * (*vec2++);
  }
}
 
/* Returns 2-norm of a double n-vector over range. */
double norm(double *vec, int end) {
  return sqrt(vectors_inner_product(end + 1, vec, vec));
}
 
void orthog1(int n, double *vec /* vector to be orthogonalized against 1 */
) {
  int i;
  double *pntr;
  double sum;
 
  sum = 0.0;
  pntr = vec;
  for (i = n; i; i--) {
    sum += *pntr++;
  }
  sum /= n;
  pntr = vec;
  for (i = n; i; i--) {
    *pntr++ -= sum;
  }
}
 
#define RANGE 500
 
void init_vec_orth1(int n, double *vec) {
  /* randomly generate a vector orthogonal to 1 (i.e., with mean 0) */
  int i;
 
  for (i = 0; i < n; i++)
    vec[i] = rand() % RANGE;
 
  orthog1(n, vec);
}
 
void right_mult_with_vector(const vtx_data *matrix, int n, const double *vector,
                            double *restrict result) {
  for (int i = 0; i < n; i++) {
    double res = 0;
    for (size_t j = 0; j < matrix[i].nedges; j++)
      res += matrix[i].ewgts[j] * vector[matrix[i].edges[j]];
    result[i] = res;
  }
}
 
void vectors_subtraction(int n, double *vector1, double *vector2,
                         double *result) {
  int i;
  for (i = 0; i < n; i++) {
    result[i] = vector1[i] - vector2[i];
  }
}
 
void vectors_addition(int n, double *vector1, double *vector2, double *result) {
  int i;
  for (i = 0; i < n; i++) {
    result[i] = vector1[i] + vector2[i];
  }
}
 
void vectors_scalar_mult(int n, const double *vector, double alpha,
                         double *result) {
  int i;
  for (i = 0; i < n; i++) {
    result[i] = vector[i] * alpha;
  }
}
 
void copy_vector(int n, const double *restrict source, double *restrict dest) {
  int i;
  for (i = 0; i < n; i++)
    dest[i] = source[i];
}
 
double vectors_inner_product(int n, const double *vector1,
                             const double *vector2) {
  int i;
  double result = 0;
  for (i = 0; i < n; i++) {
    result += vector1[i] * vector2[i];
  }
 
  return result;
}
 
double max_abs(int n, double *vector) {
  double max_val = -1e50;
  int i;
  for (i = 0; i < n; i++)
    max_val = fmax(max_val, fabs(vector[i]));
 
  return max_val;
}
 
void right_mult_with_vector_transpose(double **matrix, int dim1, int dim2,
                                      double *vector, double *result) {
  // matrix is dim2 × dim1, vector has dim2 components,
  // result = matrixᵀ × vector
  int i, j;
 
  double res;
  for (i = 0; i < dim1; i++) {
    res = 0;
    for (j = 0; j < dim2; j++)
      res += matrix[j][i] * vector[j];
    result[i] = res;
  }
}
 
void right_mult_with_vector_d(double *const *matrix, int dim1, int dim2,
                              const double *vector, double *restrict result) {
  // matrix is dim1 × dim2, vector has dim2 components,
  // result = matrix × vector
  for (int i = 0; i < dim1; i++) {
    double res = 0;
    for (int j = 0; j < dim2; j++)
      res += matrix[i][j] * vector[j];
    result[i] = res;
  }
}
 
/*****************************
** Single precision (float) **
** version                  **
*****************************/
 
void orthog1f(int n, float *vec) {
  int i;
  float *pntr;
  float sum;
 
  sum = 0.0;
  pntr = vec;
  for (i = n; i; i--) {
    sum += *pntr++;
  }
  sum /= n;
  pntr = vec;
  for (i = n; i; i--) {
    *pntr++ -= sum;
  }
}
 
void right_mult_with_vector_ff(const float *packed_matrix, int n,
                               const float *vector, float *restrict result) {
  /* packed matrix is the upper-triangular part of a symmetric matrix arranged
   * in a vector row-wise */
  for (int i = 0; i < n; i++) {
    result[i] = 0;
  }
  for (int index = 0, i = 0; i < n; i++) {
    const float vector_i = vector[i];
    /* deal with main diag */
    float res = packed_matrix[index++] * vector_i;
    /* deal with off diag */
    for (int j = i + 1; j < n; j++, index++) {
      res += packed_matrix[index] * vector[j];
      result[j] += packed_matrix[index] * vector_i;
    }
    result[i] += res;
  }
}
 
void vectors_subtractionf(int n, float *vector1, float *vector2,
                          float *result) {
  int i;
  for (i = 0; i < n; i++) {
    result[i] = vector1[i] - vector2[i];
  }
}
 
void vectors_additionf(int n, float *vector1, float *vector2, float *result) {
  int i;
  for (i = 0; i < n; i++) {
    result[i] = vector1[i] + vector2[i];
  }
}
 
void vectors_mult_additionf(int n, float *vector1, float alpha,
                            float *vector2) {
  int i;
  for (i = 0; i < n; i++) {
    vector1[i] = vector1[i] + alpha * vector2[i];
  }
}
 
void copy_vectorf(int n, float *source, float *dest) {
  int i;
  for (i = 0; i < n; i++)
    dest[i] = source[i];
}
 
double vectors_inner_productf(int n, float *vector1, float *vector2) {
  int i;
  double result = 0;
  for (i = 0; i < n; i++) {
    result += vector1[i] * vector2[i];
  }
 
  return result;
}
 
void set_vector_valf(int n, float val, float *result) {
  int i;
  for (i = 0; i < n; i++)
    result[i] = val;
}
 
double max_absf(int n, float *vector) {
  int i;
  float max_val = -1e30f;
  for (i = 0; i < n; i++)
    max_val = fmaxf(max_val, fabsf(vector[i]));
 
  return max_val;
}
 
void square_vec(int n, float *vec) {
  int i;
  for (i = 0; i < n; i++) {
    vec[i] *= vec[i];
  }
}
 
void invert_vec(int n, float *vec) {
  int i;
  for (i = 0; i < n; i++) {
    if (vec[i] != 0.0) {
      vec[i] = 1.0f / vec[i];
    }
  }
}
 
void sqrt_vecf(int n, float *source, float *target) {
  int i;
  for (i = 0; i < n; i++) {
    if (source[i] >= 0.0) {
      target[i] = sqrtf(source[i]);
    }
  }
}
 
void invert_sqrt_vec(int n, float *vec) {
  int i;
  for (i = 0; i < n; i++) {
    if (vec[i] > 0.0) {
      vec[i] = 1.0f / sqrtf(vec[i]);
    }
  }
}