/* Copyright © 2002, University of Tennessee Research Foundation. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the University of Tennessee nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ #include #include #include #include #include #include #include #include "svdlib.h" enum algorithms{LAS2}; /*********************************************************************** * * * timer() * * Returns elapsed cpu time (float, in seconds) * * * ***********************************************************************/ float timer(void) { long elapsed_time; struct rusage mytime; getrusage(RUSAGE_SELF,&mytime); /* convert elapsed time to milliseconds */ elapsed_time = (mytime.ru_utime.tv_sec * 1000 + mytime.ru_utime.tv_usec / 1000); /* return elapsed time in seconds */ return((float)elapsed_time/1000.); } static long imin(long a, long b) {return (a < b) ? a : b;} static void debug(char *fmt, ...) { va_list args; va_start(args, fmt); vfprintf(stderr, fmt, args); va_end(args); } static void fatalError(char *fmt, ...) { va_list args; va_start(args, fmt); fprintf(stderr, "ERROR: "); vfprintf(stderr, fmt, args); fprintf(stderr, "\a\n"); va_end(args); exit(1); } void printUsage(char *progname) { debug("SVD Version %s\n" "written by Douglas Rohde based on code adapted from SVDPACKC\n\n", SVDVersion); debug("usage: %s [options] matrix_file\n", progname); debug(" -a algorithm Sets the algorithm to use. They include:\n" " las2 (default)\n" " -c infile outfile\n" " Convert a matrix file to a new format (using -r and -w)\n" " Then exit immediately\n" " -d dimensions Desired SVD triples (default is all)\n" " -e bound Minimum magnitude of wanted eigenvalues (1e-30)\n" " -k kappa Accuracy parameter for las2 (1e-6)\n" " -i iterations Algorithm iterations\n" " -o file_root Root of files in which to store resulting U,S,V\n" " -r format Input matrix file format\n" " sth SVDPACK Harwell-Boeing text format\n" " st Sparse text (default)\n" " dt Dense text\n" " sb Sparse binary\n" " db Dense binary\n" " -v verbosity Default 1. 0 for no feedback, 2 for more\n" " -w format Output matrix file format (see -r for formats)\n" " (default is dense text)\n"); exit(1); } int main(int argc, char *argv[]) { extern char *optarg; extern int optind; int opt; SVDRec R = NULL; SMat A = NULL; char transpose = FALSE; int readFormat = SVD_F_ST; int writeFormat = SVD_F_DT; int algorithm = LAS2; int iterations = 0; int dimensions = 0; char *vectorFile = NULL; double las2end[2] = {-1.0e-30, 1.0e-30}; double kappa = 1e-6; double exetime; while ((opt = getopt(argc, argv, "a:c:d:e:hk:i:o:r:tv:w:")) != -1) { switch (opt) { case 'a': if (!strcasecmp(optarg, "las2")) algorithm = LAS2; else fatalError("unknown algorithm: %s", optarg); break; case 'c': if (optind != argc - 1) printUsage(argv[0]); if (SVDVerbosity > 0) printf("Converting %s to %s\n", optarg, argv[optind]); if (SVD_IS_SPARSE(readFormat) && SVD_IS_SPARSE(writeFormat)) { SMat S = svdLoadSparseMatrix(optarg, readFormat); if (!S) fatalError("failed to read sparse matrix"); if (transpose) { if (SVDVerbosity > 0) printf(" Transposing the matrix...\n"); S = svdTransposeS(S); /* This leaks memory. */ } svdWriteSparseMatrix(S, argv[optind], writeFormat); } else { DMat D = svdLoadDenseMatrix(optarg, readFormat); if (!D) fatalError("failed to read dense matrix"); if (transpose) { if (SVDVerbosity > 0) printf(" Transposing the matrix...\n"); D = svdTransposeD(D); /* This leaks memory. */ } svdWriteDenseMatrix(D, argv[optind], writeFormat); } exit(0); break; case 'd': dimensions = atoi(optarg); if (dimensions < 0) fatalError("dimensions must be non-negative"); break; case 'e': las2end[1] = atof(optarg); las2end[0] = -las2end[1]; break; case 'h': printUsage(argv[0]); break; case 'k': kappa = atof(optarg); break; case 'i': iterations = atoi(optarg); break; case 'o': vectorFile = optarg; break; case 'r': if (!strcasecmp(optarg, "sth")) { readFormat = SVD_F_STH; } else if (!strcasecmp(optarg, "st")) { readFormat = SVD_F_ST; } else if (!strcasecmp(optarg, "dt")) { readFormat = SVD_F_DT; } else if (!strcasecmp(optarg, "sb")) { readFormat = SVD_F_SB; } else if (!strcasecmp(optarg, "db")) { readFormat = SVD_F_DB; } else fatalError("bad file format: %s", optarg); break; case 't': transpose = TRUE; break; case 'v': SVDVerbosity = atoi(optarg); /*if (SVDVerbosity) printf("Verbosity = %ld\n", SVDVerbosity);*/ break; case 'w': if (!strcasecmp(optarg, "sth")) { writeFormat = SVD_F_STH; } else if (!strcasecmp(optarg, "st")) { writeFormat = SVD_F_ST; } else if (!strcasecmp(optarg, "dt")) { writeFormat = SVD_F_DT; } else if (!strcasecmp(optarg, "sb")) { writeFormat = SVD_F_SB; } else if (!strcasecmp(optarg, "db")) { writeFormat = SVD_F_DB; } else fatalError("bad file format: %s", optarg); break; default: printUsage(argv[0]); } } if (optind != argc - 1) printUsage(argv[0]); if (SVDVerbosity > 0) printf("Loading the matrix...\n"); A = svdLoadSparseMatrix(argv[optind], readFormat); if (!A) fatalError("failed to read sparse matrix. Did you specify the correct file type with the -r argument?"); if (transpose) { if (SVDVerbosity > 0) printf(" Transposing the matrix...\n"); SMat T = A; A = svdTransposeS(A); svdFreeSMat(T); } if (dimensions <= 0) dimensions = imin(A->rows, A->cols); exetime = timer(); if (SVDVerbosity > 0) printf("Computing the SVD...\n"); if (algorithm == LAS2) { if (!(R = svdLAS2(A, dimensions, iterations, las2end, kappa))) fatalError("error in svdLAS2"); } else { fatalError("unknown algorithm"); } exetime = timer() - exetime; if (SVDVerbosity > 0) { printf("\nELAPSED CPU TIME = %6g sec.\n", exetime); printf("MULTIPLICATIONS BY A = %6ld\n", (SVDCount[SVD_MXV] - R->d) / 2 + R->d); printf("MULTIPLICATIONS BY A^T = %6ld\n", (SVDCount[SVD_MXV] - R->d) / 2); } if (vectorFile) { char filename[128]; sprintf(filename, "%s-Ut", vectorFile); svdWriteDenseMatrix(R->Ut, filename, writeFormat); sprintf(filename, "%s-S", vectorFile); svdWriteDenseArray(R->S, R->d, filename, FALSE); sprintf(filename, "%s-Vt", vectorFile); svdWriteDenseMatrix(R->Vt, filename, writeFormat); } return 0; }