001    // --- BEGIN LICENSE BLOCK ---
002    /* 
003     * Copyright (c) 2009, Mikio L. Braun
004     * All rights reserved.
005     * 
006     * Redistribution and use in source and binary forms, with or without
007     * modification, are permitted provided that the following conditions are
008     * met:
009     * 
010     *     * Redistributions of source code must retain the above copyright
011     *       notice, this list of conditions and the following disclaimer.
012     * 
013     *     * Redistributions in binary form must reproduce the above
014     *       copyright notice, this list of conditions and the following
015     *       disclaimer in the documentation and/or other materials provided
016     *       with the distribution.
017     * 
018     *     * Neither the name of the Technische Universit??t Berlin nor the
019     *       names of its contributors may be used to endorse or promote
020     *       products derived from this software without specific prior
021     *       written permission.
022     * 
023     * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
024     * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
025     * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
026     * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
027     * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
028     * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
029     * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
030     * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
031     * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
032     * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
033     * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
034     */
035    // --- END LICENSE BLOCK ---
036    
037    /*
038     * To change this template, choose Tools | Templates
039     * and open the template in the editor.
040     */
041    package org.jblas.la;
042    
043    import org.jblas.la.exceptions.LapackException;
044    
045    /**
046     * <p>Implementation of some Blas functions, mostly those which require linear runtime
047     * in the number of matrix elements. Because of the copying overhead when passing
048     * primitive arrays to native code, it doesn't make sense for these functions
049     * to be implemented in native code. The Java code is about as fast.</p>
050     * 
051     * <p>The same conventions were used as in the native code, that is, for each array
052     * you also pass an index pointing to the starting index.</p>
053     * 
054     * <p>These methods are mostly optimized for the case where the starting index is 0
055     * and the increment is 1.</p>
056     */
057    public class JavaBlas {
058    
059        /** Exchange two vectors. */
060        public static void rswap(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
061            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
062                double z;
063                for (int i = 0; i < n; i++) {
064                    z = dx[i];
065                    dx[i] = dy[i];
066                    dy[i] = z;
067                }
068            } else {
069                double z;
070                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
071                    z = dx[xi];
072                    dx[xi] = dy[yi];
073                    dy[yi] = z;
074                }
075            }
076        }
077    
078        /** Copy dx to dy. */
079        public static void rcopy(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
080            if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
081                throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
082            }
083            if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
084                throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
085            }
086            if (incx == 1 && incy == 1) {
087                System.arraycopy(dx, dxIdx, dy, dyIdx, n);
088            } else {
089                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
090                    dy[yi] = dx[xi];
091                }
092            }
093        }
094    
095        /** Compute dy <- da * dx + dy. */
096        public static void raxpy(int n, double da, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
097            if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
098                throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
099            }
100            if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
101                throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
102            }
103            
104            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
105                if (da == 1.0) {
106                    for (int i = 0; i < n; i++) {
107                        dy[i] += dx[i];
108                    }
109                } else {
110                    for (int i = 0; i < n; i++) {
111                        dy[i] += da * dx[i];
112                    }
113                }
114            } else {
115                if (da == 1.0) {
116                    for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
117                        dy[yi] += dx[xi];
118                    }
119    
120                } else {
121                    for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
122                        dy[yi] += da * dx[xi];
123                    }
124                }
125            }
126        }
127    
128        /** Computes dz <- dx + dy */
129        public static void rzaxpy(int n, double[] dz, int dzIdx, int incz, double da, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
130            if (dxIdx == 0 && incx == 1 && dyIdx == 0 && incy == 1 && dzIdx == 0 && incz == 1) {
131                if (da == 1.0) {
132                    for (int c = 0; c < n; c++)
133                        dz[c] = dx[c] + dy[c];
134                } else {
135                    for (int c = 0; c < n; c++)
136                        dz[c] = da*dx[c] + dy[c];
137                }
138            } else {
139                if (da == 1.0) {
140                    for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
141                        dz[zi] = dx[xi] + dy[yi];
142                    }
143                } else {
144                    for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
145                        dz[zi] = da*dx[xi] + dy[yi];
146                    }
147                }
148            }
149        }
150    
151        public static void rzgxpy(int n, double[] dz, double[] dx, double[] dy) {
152            for (int c = 0; c < n; c++)
153                dz[c] = dx[c] + dy[c];       
154        }
155    
156        /** Compute scalar product between dx and dy. */
157        public static double rdot(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
158            double s = 0.0;
159            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
160                for (int i = 0; i < n; i++)
161                    s += dx[i] * dy[i];
162            }
163            else {
164                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
165                    s += dx[xi] * dy[yi];
166                }
167            }
168            return s;
169        }
170    //BEGIN
171      // The code below has been automatically generated.
172      // DO NOT EDIT!
173    
174        /** Exchange two vectors. */
175        public static void rswap(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
176            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
177                float z;
178                for (int i = 0; i < n; i++) {
179                    z = dx[i];
180                    dx[i] = dy[i];
181                    dy[i] = z;
182                }
183            } else {
184                float z;
185                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
186                    z = dx[xi];
187                    dx[xi] = dy[yi];
188                    dy[yi] = z;
189                }
190            }
191        }
192    
193        /** Copy dx to dy. */
194        public static void rcopy(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
195            if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
196                throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
197            }
198            if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
199                throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
200            }
201            if (incx == 1 && incy == 1) {
202                System.arraycopy(dx, dxIdx, dy, dyIdx, n);
203            } else {
204                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
205                    dy[yi] = dx[xi];
206                }
207            }
208        }
209    
210        /** Compute dy <- da * dx + dy. */
211        public static void raxpy(int n, float da, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
212            if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
213                throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
214            }
215            if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
216                throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
217            }
218            
219            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
220                if (da == 1.0f) {
221                    for (int i = 0; i < n; i++) {
222                        dy[i] += dx[i];
223                    }
224                } else {
225                    for (int i = 0; i < n; i++) {
226                        dy[i] += da * dx[i];
227                    }
228                }
229            } else {
230                if (da == 1.0f) {
231                    for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
232                        dy[yi] += dx[xi];
233                    }
234    
235                } else {
236                    for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
237                        dy[yi] += da * dx[xi];
238                    }
239                }
240            }
241        }
242    
243        /** Computes dz <- dx + dy */
244        public static void rzaxpy(int n, float[] dz, int dzIdx, int incz, float da, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
245            if (dxIdx == 0 && incx == 1 && dyIdx == 0 && incy == 1 && dzIdx == 0 && incz == 1) {
246                if (da == 1.0f) {
247                    for (int c = 0; c < n; c++)
248                        dz[c] = dx[c] + dy[c];
249                } else {
250                    for (int c = 0; c < n; c++)
251                        dz[c] = da*dx[c] + dy[c];
252                }
253            } else {
254                if (da == 1.0f) {
255                    for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
256                        dz[zi] = dx[xi] + dy[yi];
257                    }
258                } else {
259                    for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
260                        dz[zi] = da*dx[xi] + dy[yi];
261                    }
262                }
263            }
264        }
265    
266        public static void rzgxpy(int n, float[] dz, float[] dx, float[] dy) {
267            for (int c = 0; c < n; c++)
268                dz[c] = dx[c] + dy[c];       
269        }
270    
271        /** Compute scalar product between dx and dy. */
272        public static float rdot(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
273            float s = 0.0f;
274            if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
275                for (int i = 0; i < n; i++)
276                    s += dx[i] * dy[i];
277            }
278            else {
279                for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
280                    s += dx[xi] * dy[yi];
281                }
282            }
283            return s;
284        }
285    //END
286    }