/************************************************************** * Function copies a vector x to a vector y * Based on Fortran-77 routine * from Linpack by J. Dongarra * **************************************************************/ static void svd_dcopy(int n, double[] dx, int incx, double[] dy, int incy) { svd_dcopy(n, dx, 0, incx, dy, 0, incy); }
static DMat svdTransposeD(DMat D) { int r, c; DMat N = new DMat(D.cols, D.rows); for (r = 0; r < D.rows; r++) for (c = 0; c < D.cols; c++) N.value[c][r] = D.value[r][c]; return N; }
enum storeVals {STORQ, RETRQ, STORP, RETRP};
static SMat svdTransposeS(SMat S) { int r, c, i, j; SMat N = new SMat(S.cols, S.rows, S.vals); /* Count number nz in each row. */ for (i = 0; i < S.vals; i++) N.pointr[S.rowind[i]]++; /* Fill each cell with the starting point of the previous row. */ N.pointr[S.rows] = S.vals - N.pointr[S.rows - 1]; for (r = S.rows - 1; r > 0; r--) N.pointr[r] = N.pointr[r + 1] - N.pointr[r - 1]; N.pointr[0] = 0; /* Assign the new columns and values. */ for (c = 0, i = 0; c < S.cols; c++) { for (; i < S.pointr[c + 1]; i++) { r = S.rowind[i]; j = N.pointr[r + 1]++; N.rowind[j] = c; N.value[j] = S.value[i]; } } return N; }
static void svd_error(String fmt, Object... args) { svd_beep(); System.err.print("ERROR: "); System.err.printf(fmt, args); System.err.println(); }
static void svd_fatalError(String fmt, Object... args) { svd_error(fmt, args); System.exit(1); }
/***************************************************************** * Function finds the index of element having max. absolute value* based on * FORTRAN 77 routine from Linpack by J. Dongarra * *****************************************************************/ static int svd_idamax(int n, double[] dx, int ix0, int incx) { int ix,imax; double dmax; if (n < 1) return -1; if (n == 1) return 0; if (incx == 0) return -1; ix = (incx < 0) ? ix0 + ((-n+1) * incx) : ix0; imax = ix; dmax = fabs(dx[ix]); for (int i=1; i < n; i++) { ix += incx; double dtemp = fabs(dx[ix]); if (dtemp > dmax) { dmax = dtemp; imax = ix; } } return imax; }
/********************************************************************* * Function sorts array1 and array2 into increasing order for array1 * *********************************************************************/ static void svd_dsort2(int igap, int n, double[] array1, double[] array2) { double temp; int i, j, index; if (0 == igap) return; else { for (i = igap; i < n; i++) { j = i - igap; index = i; while (j >= 0 && array1[j] > array1[index]) { temp = array1[j]; array1[j] = array1[index]; array1[index] = temp; temp = array2[j]; array2[j] = array2[index]; array2[index] = temp; j -= igap; index = j + igap; } } } svd_dsort2(igap/2,n,array1,array2); }
/*********************************************************************** * * * write_header() * * Function writes out header of output file containing ritz values * * * ***********************************************************************/ void write_header(long iterations, long dimensions, double endl, double endr, boolean b, double kappa, long nrow, long ncol, long vals) { printf("SOLVING THE [A^TA] EIGENPROBLEM\n"); printf("NO. OF ROWS = %6d\n", nrow); printf("NO. OF COLUMNS = %6d\n", ncol); printf("NO. OF NON-ZERO VALUES = %6d\n", vals); printf("MATRIX DENSITY = %6.2f%%\n", ((float) vals / nrow) * 100 / ncol); /* printf("ORDER OF MATRIX A = %5ld\n", n); */ printf("MAX. NO. OF LANCZOS STEPS = %6d\n", iterations); printf("MAX. NO. OF EIGENPAIRS = %6d\n", dimensions); printf("LEFT END OF THE INTERVAL = %9.2E\n", endl); printf("RIGHT END OF THE INTERVAL = %9.2E\n", endr); printf("KAPPA = %9.2E\n", kappa); /* printf("WANT S-VECTORS? [T/F] = %c\n", (vectors) ? 'T' : 'F'); */ printf("\n"); return; }
enum storeVals {STORQ, RETRQ, STORP, RETRP};
static SMat svdConvertDtoS(DMat D) { SMat S; int i, j, n; // n = number of non-zero elements for (i = 0, n = 0; i < D.rows; i++) { for (j = 0; j < D.cols; j++) { if (D.value[i][j] != 0) n++; } } S = new SMat(D.rows, D.cols, n); for (j = 0, n = 0; j < D.cols; j++) { S.pointr[j] = n; for (i = 0; i < D.rows; i++) if (D.value[i][j] != 0) { S.rowind[n] = i; S.value[n] = D.value[i][j]; n++; } } S.pointr[S.cols] = S.vals; return S; }
/************************************************************** * Function copies a vector x to a vector y * Based on Fortran-77 routine * from Linpack by J. Dongarra * **************************************************************/ static void svd_dcopy(int n, double[] dx, int incx, double[] dy, int incy) { svd_dcopy(n, dx, 0, incx, dy, 0, incy); }
static DMat svdTransposeD(DMat D) { int r, c; DMat N = new DMat(D.cols, D.rows); for (r = 0; r < D.rows; r++) for (c = 0; c < D.cols; c++) N.value[c][r] = D.value[r][c]; return N; }
static void svd_error(String fmt, Object... args) { svd_beep(); System.err.print("ERROR: "); System.err.printf(fmt, args); System.err.println(); }
static void svd_fatalError(String fmt, Object... args) { svd_error(fmt, args); System.exit(1); }
/***************************************************************** * Function finds the index of element having max. absolute value* based on * FORTRAN 77 routine from Linpack by J. Dongarra * *****************************************************************/ static int svd_idamax(int n, double[] dx, int ix0, int incx) { int ix,imax; double dmax; if (n < 1) return -1; if (n == 1) return 0; if (incx == 0) return -1; ix = (incx < 0) ? ix0 + ((-n+1) * incx) : ix0; imax = ix; dmax = fabs(dx[ix]); for (int i=1; i < n; i++) { ix += incx; double dtemp = fabs(dx[ix]); if (dtemp > dmax) { dmax = dtemp; imax = ix; } } return imax; }
public static SMat svdConvertDtoS(DMat D) { SMat S; int i, j, n; // n = number of non-zero elements for (i = 0, n = 0; i < D.rows; i++) { for (j = 0; j < D.cols; j++) { if (D.value[i][j] != 0) n++; } } S = new SMat(D.rows, D.cols, n); for (j = 0, n = 0; j < D.cols; j++) { S.pointr[j] = n; for (i = 0; i < D.rows; i++) if (D.value[i][j] != 0) { S.rowind[n] = i; S.value[n] = D.value[i][j]; n++; } } S.pointr[S.cols] = S.vals; return S; }
static SMat svdTransposeS(SMat S) { int r, c, i, j; SMat N = new SMat(S.cols, S.rows, S.vals); /* Count number nz in each row. */ for (i = 0; i < S.vals; i++) N.pointr[S.rowind[i]]++; /* Fill each cell with the starting point of the previous row. */ N.pointr[S.rows] = S.vals - N.pointr[S.rows - 1]; for (r = S.rows - 1; r > 0; r--) N.pointr[r] = N.pointr[r + 1] - N.pointr[r - 1]; N.pointr[0] = 0; /* Assign the new columns and values. */ for (c = 0, i = 0; c < S.cols; c++) { for (; i < S.pointr[c + 1]; i++) { r = S.rowind[i]; j = N.pointr[r + 1]++; N.rowind[j] = c; N.value[j] = S.value[i]; } } return N; }