sddsmatrixmult.c File Reference

Detailed Description

Multiplies matrices from SDDS files and outputs the result.

This program reads two SDDS files containing matrix data, performs matrix multiplication, and writes the resulting product matrix to an output SDDS file. It supports options for data reuse, matrix commutation, output format selection, and verbosity for debugging and diagnostics.

Usage

sddsmatrixmult [<file1>] <file2> [<output>]
               [-pipe=[input][,output]]
               [-majorOrder=row|column]
               [-commute]
               [-reuse]
               [-verbose]
               [-ascii]
               [-threads=<number>]

Options

Option	Description
-pipe	Use input and/or output pipes.
-majorOrder	Specify output in row-major or column-major order.
-commute	Swap file1 and file2 to reverse the order of multiplication.
-reuse	Reuse the last data page if a file runs out of pages.
-ascii	Output file in ASCII mode.
-verbose	Enable detailed output for diagnostics.
-threads	Number of threads for matrix multiplication.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, L. Emery, R. Soliday, H. Shang

Definition in file sddsmatrixmult.c.

#include "mdb.h"
#include "scan.h"
#include "match_string.h"
#include "matlib.h"
#include "SDDS.h"

Go to the source code of this file.

Functions
static int	threaded_matrix_mult (MATRIX *C, MATRIX *A, MATRIX *B, int threads)
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 101 of file sddsmatrixmult.c.

                                {
  SCANNED_ARG *s_arg;
  SDDS_TABLE input1Page, input2Page, outputPage;
 
  char *inputfile1, *inputfile2, *outputfile;
  char **Input1Column, **Input1DoubleColumn;
  char **Input2Column, **Input2DoubleColumn;
  char **OutputDoubleColumn;
  long Input1Rows, Input1DoubleColumns;
  long Input2Rows, Input2DoubleColumns;
  long OutputRows;
  int32_t Input1Columns, Input2Columns, OutputDoubleColumns;
 
  long i, i_arg, col, commute;
#define BUFFER_SIZE_INCREMENT 20
  MATRIX *R1, *R1Trans, *R2, *R2Trans, *R3, *R3Trans;
  long verbose;
  long ipage, ipage1, ipage2, lastPage1, lastPage2, ascii;
  unsigned long pipeFlags, majorOrderFlag;
  long tmpfile_used, noWarnings;
  long reuse;
  int threads;
  short columnMajorOrder = -1;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&s_arg, argc, argv);
  if (argc == 1)
    bomb(NULL, USAGE);
 
  inputfile1 = inputfile2 = outputfile = NULL;
  Input1DoubleColumn = Input2DoubleColumn = OutputDoubleColumn = NULL;
  Input1Rows = Input1DoubleColumns = Input2Rows = Input2DoubleColumns = OutputRows = lastPage1 = lastPage2 = 0;
  tmpfile_used = 0;
  verbose = 0;
  pipeFlags = 0;
  noWarnings = 0;
  reuse = commute = ascii = 0;
  threads = 1;
 
  for (i_arg = 1; i_arg < argc; i_arg++) {
    if (s_arg[i_arg].arg_type == OPTION) {
      switch (match_string(s_arg[i_arg].list[0], commandline_option, N_OPTIONS, UNIQUE_MATCH)) {
      case CLO_MAJOR_ORDER:
        majorOrderFlag = 0;
        s_arg[i_arg].n_items--;
        if (s_arg[i_arg].n_items > 0 &&
            (!scanItemList(&majorOrderFlag, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
                           "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
                           "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL))) {
          SDDS_Bomb("Invalid -majorOrder syntax/values");
        }
        if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
          columnMajorOrder = 1;
        else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
          columnMajorOrder = 0;
        break;
      case CLO_PIPE:
        if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
          SDDS_Bomb("Invalid -pipe syntax");
        break;
      case CLO_VERBOSE:
        verbose = 1;
        break;
      case CLO_ASCII:
        ascii = 1;
        break;
      case CLO_REUSE:
        reuse = 1;
        break;
      case CLO_COMMUTE:
        commute = 1;
        break;
      case CLO_THREADS:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%d", &threads) != 1 || threads < 1)
          SDDS_Bomb("Invalid -threads syntax");
        break;
      default:
        bomb("Unrecognized option given", USAGE);
      }
    } else {
      if (!inputfile1)
        inputfile1 = s_arg[i_arg].list[0];
      else if (!inputfile2)
        inputfile2 = s_arg[i_arg].list[0];
      else if (!outputfile)
        outputfile = s_arg[i_arg].list[0];
      else
        bomb("Too many filenames given", USAGE);
    }
  }
 
  if (pipeFlags & USE_STDIN && inputfile1) {
    if (outputfile)
      SDDS_Bomb("Too many filenames (sddsxref)");
    outputfile = inputfile2;
    inputfile2 = inputfile1;
    inputfile1 = NULL;
  }
 
  processFilenames("sddsmatrixmult", &inputfile1, &outputfile, pipeFlags, noWarnings, &tmpfile_used);
  if (!inputfile2)
    SDDS_Bomb("Second input file not specified");
 
  if (commute) {
    char *ptr = inputfile1;
    inputfile1 = inputfile2;
    inputfile2 = ptr;
  }
 
  /* Initialize input files */
  if (!SDDS_InitializeInput(&input1Page, inputfile1))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  Input1Column = (char **)SDDS_GetColumnNames(&input1Page, &Input1Columns);
 
  if (!SDDS_InitializeInput(&input2Page, inputfile2))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  Input2Column = (char **)SDDS_GetColumnNames(&input2Page, &Input2Columns);
 
  if (!SDDS_InitializeOutput(&outputPage, ascii ? SDDS_ASCII : SDDS_BINARY, 1, "Matrix product", "Matrix product", outputfile))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  if (columnMajorOrder != -1)
    outputPage.layout.data_mode.column_major = columnMajorOrder;
  else
    outputPage.layout.data_mode.column_major = input1Page.layout.data_mode.column_major;
 
  /* Process data pages */
  while ((ipage1 = SDDS_ReadTable(&input1Page)) && (ipage2 = SDDS_ReadTable(&input2Page))) {
    if ((reuse && (ipage1 < 0 && ipage2 < 0)) ||
        (!reuse && (ipage1 < 0 || ipage2 < 0)))
      break;
 
    ipage = MAX(ipage1, ipage2);
 
    /* Process first input file */
    if (ipage1 == 1) {
      Input1DoubleColumns = 0;
      Input1DoubleColumn = (char **)malloc(Input1Columns * sizeof(char *));
      /* Count the numerical columns in input1 */
      for (i = 0; i < Input1Columns; i++) {
        if (SDDS_NUMERIC_TYPE(SDDS_GetColumnType(&input1Page, i))) {
          Input1DoubleColumn[Input1DoubleColumns] = Input1Column[i];
          Input1DoubleColumns++;
        }
      }
      if (!Input1DoubleColumns) {
        if (verbose) {
          fprintf(stderr, "No numerical columns in page %ld of file %s.\n", ipage, inputfile1 ? inputfile1 : "stdin");
        }
      }
      Input1Rows = SDDS_CountRowsOfInterest(&input1Page);
      if (Input1DoubleColumns && Input1Rows)
        m_alloc(&R1, Input1DoubleColumns, Input1Rows);
      else if (!Input1Rows) {
        if (verbose)
          fprintf(stderr, "No rows in page %ld of file %s.\n", ipage, inputfile1 ? inputfile1 : "stdin");
      }
    }
 
    if (ipage1 > 0) {
      if (Input1Rows != SDDS_CountRowsOfInterest(&input1Page))
        fprintf(stderr, "Number of rows in page %ld of file %s changed.\n", ipage, inputfile1 ? inputfile1 : "stdin");
      for (col = 0; col < Input1DoubleColumns; col++) {
        if (!(R1->a[col] = (double *)SDDS_GetColumnInDoubles(&input1Page, Input1DoubleColumn[col])))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      }
      lastPage1 = ipage1;
      if (verbose)
        fprintf(stderr, "Using page %ld of file %s.\n", lastPage1, inputfile1 ? inputfile1 : "stdin");
    } else if (ipage1 < 0) {
      if (verbose)
        fprintf(stderr, "Reusing page %ld of file %s.\n", lastPage1, inputfile1 ? inputfile1 : "stdin");
    }
 
    if (verbose && Input1DoubleColumns && Input1Rows) {
      m_alloc(&R1Trans, Input1Rows, Input1DoubleColumns);
      m_trans(R1Trans, R1);
      m_show(R1Trans, "%9.6le ", "Input matrix 1:\n", stderr);
      m_free(&R1Trans);
    }
 
    /* Process second input file */
    if (ipage2 == 1) {
      Input2DoubleColumns = 0;
      Input2DoubleColumn = (char **)malloc(Input2Columns * sizeof(char *));
      /* Count the numerical columns in input2 */
      for (i = 0; i < Input2Columns; i++) {
        if (SDDS_NUMERIC_TYPE(SDDS_GetColumnType(&input2Page, i))) {
          Input2DoubleColumn[Input2DoubleColumns] = Input2Column[i];
          Input2DoubleColumns++;
        }
      }
      if (!Input2DoubleColumns) {
        if (verbose) {
          fprintf(stderr, "No numerical columns in page %ld of file %s.\n", ipage, inputfile2);
        }
      }
      Input2Rows = SDDS_CountRowsOfInterest(&input2Page);
      if (Input2DoubleColumns && Input2Rows)
        m_alloc(&R2, Input2DoubleColumns, Input2Rows);
      else if (!Input2Rows) {
        if (verbose)
          fprintf(stderr, "No rows in page %ld of file %s.\n", ipage, inputfile2);
      }
    }
 
    if (ipage2 > 0) {
      if (Input2Rows != SDDS_CountRowsOfInterest(&input2Page))
        fprintf(stderr, "Number of rows in page %ld of file %s changed.\n", ipage, inputfile2);
      for (col = 0; col < Input2DoubleColumns; col++) {
        if (!(R2->a[col] = (double *)SDDS_GetColumnInDoubles(&input2Page, Input2DoubleColumn[col])))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      }
      lastPage2 = ipage2;
      if (verbose)
        fprintf(stderr, "Using page %ld of file %s.\n", lastPage2, inputfile2 ? inputfile2 : "stdin");
    } else if (ipage2 < 0) {
      if (verbose)
        fprintf(stderr, "Reusing page %ld of file %s.\n", lastPage2, inputfile2 ? inputfile2 : "stdin");
    }
 
    if (verbose && Input2DoubleColumns && Input2Rows) {
      m_alloc(&R2Trans, Input2Rows, Input2DoubleColumns);
      m_trans(R2Trans, R2);
      m_show(R2Trans, "%9.6le ", "Input matrix 2:\n", stderr);
      m_free(&R2Trans);
    }
 
    /* Define output table */
    if (ipage == 1) {
      OutputRows = Input1Rows;
      OutputDoubleColumns = Input2DoubleColumns;
      if (Input1DoubleColumns != Input2Rows) {
        fprintf(stderr, "Error: Dimension mismatch in files.\n");
        fprintf(stderr, "Right-hand matrix (%s) is %ldx%ld.\n",
                inputfile2 ? inputfile2 : "stdin", Input2Rows, Input2DoubleColumns);
        fprintf(stderr, "Left-hand matrix (%s) is %ldx%ld.\n",
                inputfile1 ? inputfile1 : "stdin", Input1Rows, Input1DoubleColumns);
        exit(EXIT_FAILURE);
      }
    }
 
    /* Perform matrix multiplication */
    if (OutputRows && OutputDoubleColumns) {
      if (ipage == 1)
        m_alloc(&R3, OutputDoubleColumns, OutputRows);
      if (verbose)
        fprintf(stderr, "Multiplying %d x %d matrix by %d x %d matrix\n", R2->m, R2->n, R1->m, R1->n);
      if (!threaded_matrix_mult(R3, R2, R1, threads))
        SDDS_Bomb("Matrix multiplication failed");
      if (verbose) {
        m_alloc(&R3Trans, OutputRows, OutputDoubleColumns);
        m_trans(R3Trans, R3);
        m_show(R3Trans, "%9.6le ", "Output matrix:\n", stderr);
        m_free(&R3Trans);
      }
    } else {
      if (verbose)
        fprintf(stderr, "Output file will either have no columns or no rows in page %ld.\n", ipage);
    }
 
    if (ipage == 1) {
      for (i = 0; i < Input2DoubleColumns; i++) {
        if (SDDS_TransferColumnDefinition(&outputPage, &input2Page, Input2DoubleColumn[i], NULL) < 0)
          SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
      }
      OutputDoubleColumn = (char **)SDDS_GetColumnNames(&outputPage, &OutputDoubleColumns);
      if (!SDDS_WriteLayout(&outputPage))
        SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
    }
 
    if (!SDDS_StartTable(&outputPage, OutputRows))
      SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
 
    /* Assign values to output table */
    if (OutputRows && OutputDoubleColumns) {
      for (i = 0; i < OutputDoubleColumns; i++) { /* i is the column index */
        if (!SDDS_SetColumnFromDoubles(&outputPage, SDDS_SET_BY_NAME | SDDS_PASS_BY_REFERENCE,
                                       R3->a[i], OutputRows, OutputDoubleColumn[i]))
          SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
      }
    }
 
    if (!SDDS_WriteTable(&outputPage))
      SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
  }
 
  /* Terminate SDDS tables */
  if (!SDDS_Terminate(&input1Page) || !SDDS_Terminate(&input2Page) || !SDDS_Terminate(&outputPage))
    SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
 
  if (tmpfile_used && !replaceFileAndBackUp(inputfile1, outputfile))
    return EXIT_FAILURE;
 
  return EXIT_SUCCESS;
}

threaded_matrix_mult()

static int threaded_matrix_mult ( MATRIX * C, MATRIX * A, MATRIX * B, int threads )

static

Definition at line 401 of file sddsmatrixmult.c.

                                                                              {
  long i;
  long n, m, p;
 
  if ((m = A->m) != B->n || (n = A->n) != C->n || (p = B->m) != C->m)
    return 0;
  if (threads <= 1)
    return m_mult(C, A, B);
 
#pragma omp parallel for if (threads > 1) num_threads(threads)
  for (i = 0; i < n; i++) {
    long j, k;
    double *a_i = A->a[i];
    double *c_i = C->a[i];
    for (j = 0; j < p; j++) {
      double sum = 0;
      for (k = 0; k < m; k++)
        sum += a_i[k] * B->a[k][j];
      c_i[j] = sum;
    }
  }
  return 1;
}