001 package numal.highprecision; 002 003 import numal.highprecision.*; 004 005 /** A library of basic operations between matrices and vectrors. 006 * These methods make use of lower and upper indices; which allows 007 * vectors to be 0-index or 1-index (i.e. a[0...n-1] or a[1...n]) 008 */ 009 public class Basic extends Object { 010 011 static final int BASE = 100; 012 013 /** Compute a constant multiple (<i>x</i) of part of a vector <i>b</i>. 014 * Method originally from class numal.LinearAlgebra - Lau (2004) 015 * section 1.2.A Real vector and matrix - Duplication p.6. 016 * 017 * @param l lower row bound index 018 * @param u upper row bound index 019 * @param shift index shifting parameter 020 * @param a rectangular output matrix 021 * @param b rectangular input matrix 022 */ 023 public static void dupvec(int l, int u, int shift, double a[], double b[]) { 024 for (; l<=u; l++) a[l]=b[l+shift]; 025 } 026 027 /** Replace a column sequence of elements of a rectangular matrix <i>a</i> 028 * by a constant multiple (<i>x</i) of a column sequence of elements of a rectangular matrix <i>b</i>. 029 * Method originally from class numal.LinearAlgebra - Lau (2004) 030 * section 1.3.C Real vector and matrix - Multiplication p.9. 031 * 032 * @param l lower row bound 033 * @param u upper row bound 034 * @param i column index of <i>a</i> 035 * @param j column index of <i>b</i> 036 * @param a rectangular matrix 037 * @param b rectangular matrix 038 * @param x multiplication factor 039 */ 040 public static void mulcol(int l, int u, int i, int j, double a[][], double b[][], double x) { 041 for (; l<=u; l++) a[l][i]=b[l][j]*x; 042 } 043 044 /** Compute the inner product of part of a vector <i>a</i> and part of a vector <i>b</i> 045 * A subset of the inner product 046 * can be computed by choosing suitable lower and upper bounds. 047 * Method originally from class numal.Basic - Lau (2004) 048 * section 1.4.A Real vector vector products p.10. 049 * 050 * @param l lower bound 051 * @param u upper bound 052 * @param shift index-shifting parameter of vector <i>b</i> 053 * @param a vector 054 * @param b vector 055 * @return the inner product of vector <i>a</i> and vector <i>b</i> 056 */ 057 public static double vecvec(int l, int u, int shift, double a[], double b[]) { 058 int k; 059 double s; 060 061 s=0.0; 062 for (k=l; k<=u; k++) s += a[k]*b[k+shift]; 063 return (s); 064 } 065 066 /** Compute the inner product of part of a row of a rectangular matrix <i>a</i> 067 * and the corresponding part of a vector <i>b</i>. A subset of the inner product 068 * can be computed by choosing suitable lower and upper column bounds. 069 * Method originally from class numal.Basic - Lau (2004) 070 * section 1.4.B Real Vector Vector products p.11. 071 * 072 * @param l lower bound 073 * @param u upper bound 074 * @param i row index of <i>a</i> 075 * @param a rectangular matrix 076 * @param b vector 077 * @return the inner product of row <i>i</i> of matrix <i>a</i> with vector <i>b</i> 078 */ 079 public static double matvec(int l, int u, int i, double a[][], double b[]) { 080 int k; 081 double s; 082 083 s=0.0; 084 for (k=l; k<=u; k++) s += a[i][k]*b[k]; 085 return (s); 086 } 087 088 /** 089 * Method originally from class numal.Basic - Lau (2004) 090 * section 1.4.B Real Vector Vector products p.12. 091 * @param l lower bound 092 * @param u upper bound 093 * @param i column index of <i>a</i> 094 * @param a rectangular matrix with column [l-u][j] defined 095 * @param b vector with entries [l-u] defined 096 */ 097 public static double tamvec(int l, int u, int i, double a[][], double b[]) { 098 int k; 099 double s; 100 101 s=0.0; 102 for (k=l; k<=u; k++) s += a[k][i]*b[k]; 103 return (s); 104 } 105 106 /** Get the inner product of part of a row of a rectangular matrix <i>a</i> 107 * and the corresponding part of a column of a rectangular matrix <i>b</i>. 108 * 109 * @param l lower row bound 110 * @param u upper row bound 111 * @param i row index of <i>a</i> 112 * @param j column index of <i>b</i> 113 * @param a rectangular matrix with row [i][l-u] defined 114 * @param b rectangular matrix with column [l-u][j] defined 115 * 116 * @return inner product of <i>a</i> and <i>b</i> 117 */ 118 public static double matmat(int l, int u, int i, int j, double a[][], double b[][]) { 119 int k; 120 double s; 121 122 s=0.0; 123 for (k=l; k<=u; k++) s += a[i][k]*b[k][j]; 124 return (s); 125 } 126 127 /** Compute the inner product of part of a column of a rectangular matrix <i>a</i> 128 * and the corresponding part of a column of a rectangular matrix <i>b</i>. 129 * Method originally from class numal.Basic - Lau (2004) 130 * section 1.4.E Real Vector Vector products p.13. 131 * 132 * @param l lower row bound 133 * @param u upper row bound 134 * @param i column index of <i>a</i> 135 * @param j column index of <i>b</i> 136 * @param a rectangular matrix with column [l-u][i] defined 137 * @param b rectangular matrix with column [l-u][j] defined 138 */ 139 public static double tammat(int l, int u, int i, int j, double a[][], double b[][]) { 140 int k; 141 double s; 142 143 s=0.0; 144 for (k=l; k<=u; k++) s += a[k][i]*b[k][j]; 145 return (s); 146 } 147 148 /** Compute the inner product of part of a row of a rectangular matrix <i>a</i> and the 149 * corresponding part of a row of a rectangular matrix <i>b</i>. 150 * Method originally from class numal.LinearAlgebra - Lau (2004) 151 * section 1.4.F Real Vector Vector products p.13 152 * 153 * @param l lower column bound 154 * @param u upper column bound 155 * @param i row index of <i>a</i> 156 * @param j row index of <i>b</i> 157 * @param a rectangular matrix 158 * @param b rectangular matrix 159 * @return inner product of row <i>i</i> of <i>a</i> and row <i>j</i> of <i>b</i> 160 */ 161 public static double mattam(int l, int u, int i, int j, double a[][], double b[][]) { 162 int k; 163 double s; 164 165 s=0.0; 166 for (k=l; k<=u; k++) s += a[i][k]*b[j][k]; 167 return (s); 168 } 169 170 /** Adds a constant multiple (<i>x</i>) of part of a column of a rectangular matrix 171 * <i>b</i> to part of a column of a rectangular matrix <i>a</i> 172 * Method originally from class numal.Basic - Lau (2004) 173 * section 1.7.B Real vector and matrix - Elimination p.22. 174 * 175 * @param l lower row bound 176 * @param u upper row bound 177 * @param i column index of <i>a</i> 178 * @param j column index of <i>b</i> 179 * @param a rectangular matrix 180 * @param b rectangular matrix 181 * @param x multiplication factor 182 */ 183 public static void elmcol(int l, int u, int i, int j, double a[][], double b[][], double x) { 184 for (; l<=u; l++) a[l][i] += b[l][j]*x; 185 } 186 187 /** Adds a constant multiple (<i>x</i>) of part of a row of a rectangular matrix 188 * <i>b</i> to part of a row of a rectangular matrix <i>a</i> 189 * Method originally from class numal.Basic - Lau (2004) 190 * section 1.7.B Real vector and matrix - Elimination p.22. 191 * 192 * @param l lower row bound 193 * @param u upper row bound 194 * @param i column index of <i>a</i> 195 * @param j column index of <i>b</i> 196 * @param a rectangular matrix 197 * @param b rectangular matrix 198 * @param x multiplication factor 199 */ 200 public static void elmrow(int l, int u, int i, int j, double a[][], double b[][], double x) { 201 for (; l<=u; l++) a[i][l] += b[j][l]*x; 202 } 203 204 /** Rotates two columns of a rectangular matrix using two parameters 205 * <ul> 206 * <li>a[k][i] = c * a[k][i] + s * a[k][j] (where k = l to u)</li> 207 * <li>a[k][j] = c * a[k][j] - s * a[k][i]</li> 208 * </ul> 209 * Method originally from class numal.Basic - Lau (2004) 210 * section 1.9.A Real vector and matrix - Rotation p.30. 211 * 212 * @param l lower row bound 213 * @param u upper row bound 214 * @param i column index of <i>a</i> 215 * @param j column index of <i>a</i> 216 * @param a rectangular matrix 217 * @param c cosine multiplication parameter 218 * @param s sine multiplication parameter 219 */ 220 public static void rotcol(int l, int u, int i, int j, double a[][], double c, double s) { 221 double x, y; 222 223 for (; l<=u; l++) { 224 x=a[l][i]; 225 y=a[l][j]; 226 a[l][i]=x*c+y*s; 227 a[l][j]=y*c-x*s; 228 } 229 } 230 231 }