tProgram validated - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit aae8e123fc8478c2af55280b6c26695d9582cd8f
(DIR) parent 59c80bc77e0f38e44b17dc1ebe3d44fb74a5f316
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 23 Jan 2013 07:49:03 +0100
Program validated
Diffstat:
M matrixmul/c-arrofarrs.c | 85 +++++++++++++++----------------
M matrixmul/c-gsl-cblas.c | 48 ++++++++++++++++----------------
M matrixmul/c-linarr.c | 80 ++++++++++++++++----------------
M matrixmul/cpp-linvectors.cpp | 71 +++++++++++++++++--------------
M matrixmul/cpp-vectorofvectors.cpp | 73 ++++++++++++++++++-------------
M matrixmul/lua-arrofarrs.lua | 26 +++++++++++++-------------
M matrixmul/octave.m | 2 +-
M matrixmul/python-numpy.py | 4 ++--
8 files changed, 202 insertions(+), 187 deletions(-)
---
(DIR) diff --git a/matrixmul/c-arrofarrs.c b/matrixmul/c-arrofarrs.c
t@@ -1,64 +1,59 @@
#include <stdio.h>
#include <stdlib.h>
-void matrixMult(double** A, double** B, double** C, int N)
+void matrixMult(double** A, double** B, double** C, unsigned int N)
{
- int i, j, k;
- double sum;
+ unsigned int i, j, k;
+ double sum;
#pragma omp parallel for private (j,k,sum) shared(A,B,C,N) default(none)
- for (i = 0; i<N; i++) {
- for (j = 0; j<N; j++) {
- sum = 0.0f;
- for (k = 0; k<N; k++) {
- sum += A[i][k] * B[k][j];
- }
- C[i][j] = sum;
+ for (i = 0; i<N; i++) {
+ for (j = 0; j<N; j++) {
+ sum = 0.0;
+ for (k = 0; k<N; k++) {
+ sum += A[i][k] * B[k][j];
+ }
+ C[i][j] = sum;
+ }
}
- }
}
int main(int argc, char* argv[])
{
- unsigned int i, j, N;
- double** A;
- double** B;
- double** C;
-
- if (argc == 2) {
- N = atoi(argv[1]);
- } else {
- printf("Sorry, I need matrix width as command line argument\n");
- return 1;
- }
+ unsigned int i, j, N;
+ double** A;
+ double** B;
+ double** C;
+
+ if (argc == 2) {
+ N = atoi(argv[1]);
+ } else {
+ fprintf(stderr, "Sorry, I need matrix width as command line argument\n");
+ return 1;
+ }
- A = (double**) malloc(N * sizeof(double*));
- B = (double**) malloc(N * sizeof(double*));
- C = (double**) malloc(N * sizeof(double*));
+ A = (double**) malloc(N * sizeof(double*));
+ B = (double**) malloc(N * sizeof(double*));
+ C = (double**) malloc(N * sizeof(double*));
- for (i = 0; i < N; i++) {
- A[i] = (double*) malloc(N * sizeof(double));
- B[i] = (double*) malloc(N * sizeof(double));
- C[i] = (double*) malloc(N * sizeof(double));
- }
+ for (i = 0; i < N; i++) {
+ A[i] = (double*) malloc(N * sizeof(double));
+ B[i] = (double*) malloc(N * sizeof(double));
+ C[i] = (double*) malloc(N * sizeof(double));
+ }
- for (i = 0; i < N; i++) {
- for (j = 0; j < N; j++) {
- A[i][j] = 2.0;
- B[i][j] = (double) N*j + i;
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < N; j++) {
+ A[i][j] = 2.0;
+ B[i][j] = (double) N*j + i;
+ }
}
- }
- matrixMult(A, B, C, N);
+ matrixMult(A, B, C, N);
- for (i = 0; i < N; i++) {
- free(A[i]);
- free(B[i]);
- free(C[i]);
- }
- free(A);
- free(B);
- free(C);
+ free(A);
+ free(B);
+ free(C);
- return 0;
+ return 0;
}
(DIR) diff --git a/matrixmul/c-gsl-cblas.c b/matrixmul/c-gsl-cblas.c
t@@ -4,34 +4,34 @@
int main(int argc, char* argv[])
{
- unsigned int i, N;
- double* A;
- double* B;
- double* C;
+ unsigned int i, N;
+ double* A;
+ double* B;
+ double* C;
- if (argc == 2) {
- N = atoi(argv[1]);
- } else {
- printf("Sorry, I need matrix width as command line argument\n");
- return 1;
- }
+ if (argc == 2) {
+ N = atoi(argv[1]);
+ } else {
+ printf("Sorry, I need matrix width as command line argument\n");
+ return 1;
+ }
- A = (double*) malloc(N * N * sizeof(double*));
- B = (double*) malloc(N * N * sizeof(double*));
- C = (double*) malloc(N * N * sizeof(double*));
+ A = (double*) malloc(N * N * sizeof(double*));
+ B = (double*) malloc(N * N * sizeof(double*));
+ C = (double*) malloc(N * N * sizeof(double*));
- for (i = 0; i < N*N; i++) {
- A[i] = 2.0;
- B[i] = (double)i;
- }
+ for (i = 0; i < N*N; i++) {
+ A[i] = 2.0;
+ B[i] = (double)i;
+ }
- cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
- N, N, N, 1.0, A, N, B, N, 0.0, C, N);
+ cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+ N, N, N, 1.0, A, N, B, N, 0.0, C, N);
- free(A);
- free(B);
- free(C);
+ free(A);
+ free(B);
+ free(C);
- /* Exit with success */
- return 0;
+ /* Exit with success */
+ return 0;
}
(DIR) diff --git a/matrixmul/c-linarr.c b/matrixmul/c-linarr.c
t@@ -1,52 +1,52 @@
#include <stdio.h>
#include <stdlib.h>
-void matrixMult(double* A, double* B, double* C, int N)
+void matrixMult(double* A, double* B, double* C, unsigned int N)
{
- int i, j, k;
- double sum;
+ unsigned int i, j, k;
+ double sum;
#pragma omp parallel for private (j,k,sum) default(none) shared(A,B,C,N) if(N>100)
- for (i = 0; i<N; i++) {
- for (j = 0; j<N; j++) {
- sum = 0.0f;
- for (k = 0; k<N; k++) {
- sum += A[k*N+i] * B[j*N+k];
- }
- C[j*N+i] = sum;
+ for (i = 0; i<N; i++) {
+ for (j = 0; j<N; j++) {
+ sum = 0.0;
+ for (k = 0; k<N; k++) {
+ sum += A[k*N+i] * B[j*N+k];
+ }
+ C[j*N+i] = sum;
+ }
}
- }
}
int main(int argc, char* argv[])
{
- unsigned int i, N;
- double* A;
- double* B;
- double* C;
-
- if (argc == 2) {
- N = atoi(argv[1]);
- } else {
- printf("Sorry, I need matrix width as command line argument\n");
- return 1;
- }
-
- A = (double*) malloc(N * N * sizeof(double*));
- B = (double*) malloc(N * N * sizeof(double*));
- C = (double*) malloc(N * N * sizeof(double*));
-
- for (i = 0; i < N*N; i++) {
- A[i] = 2.0;
- B[i] = (double)i;
- }
-
- matrixMult(A, B, C, N);
-
- free(A);
- free(B);
- free(C);
-
- /* Exit with success */
- return 0;
+ unsigned int i, N;
+ double* A;
+ double* B;
+ double* C;
+
+ if (argc == 2) {
+ N = atoi(argv[1]);
+ } else {
+ fprintf(stderr, "Sorry, I need matrix width as command line argument\n");
+ return 1;
+ }
+
+ A = (double*) malloc(N * N * sizeof(double*));
+ B = (double*) malloc(N * N * sizeof(double*));
+ C = (double*) malloc(N * N * sizeof(double*));
+
+ for (i = 0; i < N*N; i++) {
+ A[i] = 2.0;
+ B[i] = (double)i;
+ }
+
+ matrixMult(A, B, C, N);
+
+ free(A);
+ free(B);
+ free(C);
+
+ /* Exit with success */
+ return 0;
}
(DIR) diff --git a/matrixmul/cpp-linvectors.cpp b/matrixmul/cpp-linvectors.cpp
t@@ -2,40 +2,49 @@
#include <cstdlib>
#include <vector>
+using std::vector;
+
+void matrixMult(vector<double>& A, vector<double>& B, vector<double>& C, unsigned int N)
+{
+ unsigned int i, j, k;
+ double sum;
+ for (i = 0; i<N; i++) {
+ for (j = 0; j<N; j++) {
+ sum = 0.0;
+ for (k = 0; k<N; k++) {
+ sum += A[k*N+i] * B[j*N+k];
+ }
+ C[j*N+i] = sum;
+ }
+ }
+}
+
+
int main(int argc, char* argv[])
{
- using std::cout;
- using std::vector;
-
- unsigned int N, i, j, k;
-
- if (argc == 2) {
- N = atoi(argv[1]);
- } else {
- cout << "Sorry, I need matrix width as command line argument\n";
- return 1;
- }
-
- vector<double> A(N*N);
- vector<double> B(N*N);
- vector<double> C(N*N);
-
- for (i = 0; i<N; ++i) {
- for (j = 0; j<N; ++j) {
- A[j*N+i] = 2.0;
- B[j*N+i] = (double)i;
+ using std::cout;
+
+ unsigned int N, i, j;
+
+ if (argc == 2) {
+ N = atoi(argv[1]);
+ } else {
+ std::cerr << "Sorry, I need matrix width as command line argument\n";
+ return 1;
}
- }
-
- double sum;
- for (i = 0; i < N; ++i) {
- for (j = 0; j < N; ++j) {
- sum = 0.0;
- for (k = 0; k < N; ++k)
- sum += A[k*N+i] * B[j*N+k];
- C[j*N+i] = sum;
+
+ vector<double> A(N*N);
+ vector<double> B(N*N);
+ vector<double> C(N*N);
+
+ for (i = 0; i<N; ++i) {
+ for (j = 0; j<N; ++j) {
+ A[j*N+i] = 2.0;
+ B[j*N+i] = (double) N*j + i;
+ }
}
- }
- return 0;
+ matrixMult(A, B, C, N);
+
+ return 0;
}
(DIR) diff --git a/matrixmul/cpp-vectorofvectors.cpp b/matrixmul/cpp-vectorofvectors.cpp
t@@ -2,40 +2,51 @@
#include <cstdlib>
#include <vector>
+using std::vector;
+
+void matrixMult(
+ vector< vector<double> >& A,
+ vector< vector<double> >& B,
+ vector< vector<double> >& C,
+ unsigned int N)
+{
+ unsigned int i, j, k;
+ double sum;
+ for (i = 0; i < N; ++i) {
+ for (j = 0; j < N; ++j) {
+ sum = 0.0;
+ for (k = 0; k < N; ++k)
+ sum += A[i][k] * B[k][j];
+ C[i][j] = sum;
+ }
+ }
+}
+
int main(int argc, char* argv[])
{
- using std::cout;
- using std::vector;
-
- unsigned int N, i, j, k;
-
- if (argc == 2) {
- N = atoi(argv[1]);
- } else {
- cout << "Sorry, I need matrix width as command line argument\n";
- return 1;
- }
-
- vector< vector<double> > A(N,vector<double>(N));
- vector< vector<double> > B(N,vector<double>(N));
- vector< vector<double> > C(N,vector<double>(N));
-
- for (i = 0; i<N; ++i) {
- for (j = 0; j<N; ++j) {
- A[i][j] = 2.0;
- B[i][j] = (double) N*j + i;
+ using std::cout;
+
+ unsigned int N, i, j;
+
+ if (argc == 2) {
+ N = atoi(argv[1]);
+ } else {
+ std::cerr << "Sorry, I need matrix width as command line argument\n";
+ return 1;
}
- }
-
- double sum;
- for (i = 0; i < N; ++i) {
- for (j = 0; j < N; ++j) {
- sum = 0.0;
- for (k = 0; k < N; ++k)
- sum = A[k][j] * B[i][k];
- C[i][j] = sum;
+
+ vector< vector<double> > A(N,vector<double>(N));
+ vector< vector<double> > B(N,vector<double>(N));
+ vector< vector<double> > C(N,vector<double>(N));
+
+ for (i = 0; i<N; ++i) {
+ for (j = 0; j<N; ++j) {
+ A[i][j] = 2.0;
+ B[i][j] = (double) N*j + i;
+ }
}
- }
- return 0;
+ matrixMult(A, B, C, N);
+
+ return 0;
}
(DIR) diff --git a/matrixmul/lua-arrofarrs.lua b/matrixmul/lua-arrofarrs.lua
t@@ -7,21 +7,21 @@ A = {}
B = {}
C = {}
for i=1,N do
- A[i] = {}
- B[i] = {}
- C[i] = {}
- for j=1,N do
- A[i][j] = 2
- B[i][j] = (N * j-1) + i-1
- end
+ A[i] = {}
+ B[i] = {}
+ C[i] = {}
+ for j=1,N do
+ A[i][j] = 2
+ B[i][j] = N * (j-1) + i-1
+ end
end
for i=1,N do
- for j=1,N do
- sum = 0.0
- for k=1,N do
- sum = sum + A[i][k] * B[k][j]
+ for j=1,N do
+ sum = 0.0
+ for k=1,N do
+ sum = sum + A[i][k] * B[k][j]
+ end
+ C[i][j] = sum
end
- C[i][j] = sum
- end
end
(DIR) diff --git a/matrixmul/octave.m b/matrixmul/octave.m
t@@ -1,6 +1,6 @@
#/usr/bin/octave -q -f --no-window-system
-if (nargin == 1)
+if (nargin > 0)
N = str2num(argv(){1});
else
disp("Sorry, I need the matrix size as a command line argument");
(DIR) diff --git a/matrixmul/python-numpy.py b/matrixmul/python-numpy.py
t@@ -8,6 +8,6 @@ else :
print("Sorry, I need a matrix width as input argument!")
sys.exit(1)
-A = numpy.ones((N,N))*2.0
-B = numpy.arange(0,N*N).reshape(N,N)
+A = numpy.ones((N, N))*2.0
+B = numpy.arange(0.0, N*N).reshape(N, N, order='F')
C = numpy.dot(A,B)