\(A A^{-1} = E \\ \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)
µÕ¹ÔÎó¤òµá¤á¤ë ¢Î Í¿¤¨¤é¤ì¤¿ \( a_{ij} \) ¤ËÂФ¹¤ë \( x_{ij} \) ¤òµá¤á¤ëϢΩ°ì¼¡ÊýÄø¼°¤ò²ò¤¯
\(A = LU \)
¤È¤Ê¤ë L U ¤òµá¤á¤ë\(\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \\ \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \)
\(A A^{-1} = E \\ L U A^{-1} = E \\ U A^{-1} = L^{-1} E \\ A^{-1} = U^{-1} L^{-1} E //\\ \)
\(\begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \\ = \begin{bmatrix} l_{11} & l_{11}u_{12} & l_{11}u_{13} & l_{11}u_{14} \\ l_{21} & l_{21}u_{12} + l_{22} & l_{21}u_{13} + l_{22}u_{23} & l_{21}u_{14} + l_{22}u_{24} \\ l_{31} & l_{31}u_{12} + l_{32} & l_{31}u_{13} + l_{32}u_{23} + l_{33} & l_{31}u_{14} + l_{32}u_{24} + l_{33}u_{34} \\ l_{41} & l_{41}u_{12} + l_{42} & l_{41}u_{13} + l_{42}u_{23} + l_{43} & l_{41}u_{14} + l_{42}u_{24} + l_{43}u_{34} + l_{44} \end{bmatrix} \\ = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \\ \end{bmatrix} \)
\({\color{red}{l_{11}}} = a_{11} \\ l_{11}{\color{red}{u_{12}}} = a_{12} \\ l_{11}{\color{red}{u_{13}}} = a_{13} \\ l_{11}{\color{red}{u_{14}}} = a_{14} \\ {\color{red}{l_{21}}} = a_{21} \\ l_{21}u_{12} + {\color{red}{l_{22}}} = a_{22} \\ l_{21}u_{13} + l_{22}{\color{red}{u_{23}}} = a_{23} \\ l_{21}u_{14} + l_{22}{\color{red}{u_{24}}} = a_{24} \\ {\color{red}{l_{31}}} = a_{31} \\ l_{31}u_{12} + {\color{red}{l_{32}}} = a_{32} \\ l_{31}u_{13} + l_{32}u_{23} + {\color{red}{l_{33}}} = a_{33} \\ l_{31}u_{14} + l_{32}u_{24} + l_{33}{\color{red}{u_{34}}} = a_{34} \\ {\color{red}{l_{41}}} = a_{41} \\ l_{41}u_{12} + {\color{red}{l_{42}}} = a_{42} \\ l_{41}u_{13} + l_{42}u_{23} + {\color{red}{l_{43}}} = a_{43} \\ l_{41}u_{14} + l_{42}u_{24} + l_{43}u_{34} + {\color{red}{l_{44}}} = a_{44} \\ \)
(ÀÖ»ú ¤¬¤½¤Î¼°¤òɾ²Á¤·¤Æ¤¤¤ë»þÅÀ¤Ç¤Î̤Ãοô ¢ª ¤É¤Î¼°¤Î̤Ãοô¤â£±¤Ä¤Ê¤Î¤Çµá¤á¤ë¤³¤È¤¬¤Ç¤¤ë)\(ie. \\ {\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ {\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
\({\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ {\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
\(A A^{-1} = E \\ \begin{bmatrix} {\color{red}{a_{11}}} & {\color{red}{a_{12}}} & {\color{red}{a_{13}}} & {\color{red}{a_{14}}} \\ {\color{green}{a_{21}}} & {\color{green}{a_{22}}} & {\color{green}{a_{23}}} & {\color{green}{a_{24}}} \\ {\color{blue}{a_{31}}} & {\color{blue}{a_{32}}} & {\color{blue}{a_{33}}} & {\color{blue}{a_{34}}} \\ {\color{pink}{a_{41}}} & {\color{pink}{a_{42}}} & {\color{pink}{a_{43}}} & {\color{pink}{a_{44}}} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} \)
\(A' A^{-1} = B \\ \begin{bmatrix} {\color{blue}{a_{31}}} & {\color{blue}{a_{32}}} & {\color{blue}{a_{33}}} & {\color{blue}{a_{34}}} \\ {\color{green}{a_{21}}} & {\color{green}{a_{22}}} & {\color{green}{a_{23}}} & {\color{green}{a_{24}}} \\ {\color{red}{a_{11}}} & {\color{red}{a_{12}}} & {\color{red}{a_{13}}} & {\color{red}{a_{14}}} \\ {\color{pink}{a_{41}}} & {\color{pink}{a_{42}}} & {\color{pink}{a_{43}}} & {\color{pink}{a_{44}}} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} \)
\(A A^{-1} = E \\ A' A^{-1} = B \\ L' U' A^{-1} = B \\ U' A^{-1} = L'^{-1} B \\ A^{-1} = U'^{-1} L'^{-1} B \\ \)
°ìÈ̤ΰ켡ÊýÄø¼°¤ò²ò¤¯¤È¤¤Ë¤Ï¡¢B ¤ò³Ý¤±¤ëɬÍפ¬¤¢¤ë¤±¤É¡¢µÕ¹ÔÎó¤òµá¤á¤ë¤È¤¸ÂÄê¤Ç \(Y = U'^{-1} L'^{-1}\) ¤ÎÎó¤òÆþ¤ìÂØ¤¨¤Æ¤ä¤ë¤À¤±¤Ç¤¤¤¤\(A^{-1} = U'^{-1} L'^{-1} B = Y B = \\ \begin{bmatrix} y_{11} & y_{12} & y_{13} & y_{14} \\ y_{21} & y_{22} & y_{23} & y_{24} \\ y_{31} & y_{32} & y_{33} & y_{34} \\ y_{41} & y_{42} & y_{43} & y_{44} \end{bmatrix} \begin{bmatrix} {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} = \begin{bmatrix} {\color{red}{y_{13}}} & {\color{green}{y_{12}}} & {\color{blue}{y_{11}}} & {\color{pink}{y_{14}}} \\ {\color{red}{y_{23}}} & {\color{green}{y_{22}}} & {\color{blue}{y_{21}}} & {\color{pink}{y_{24}}} \\ {\color{red}{y_{33}}} & {\color{green}{y_{32}}} & {\color{blue}{y_{31}}} & {\color{pink}{y_{34}}} \\ {\color{red}{y_{43}}} & {\color{green}{y_{42}}} & {\color{blue}{y_{41}}} & {\color{pink}{y_{44}}} \end{bmatrix} \)
\({\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ \)
\({\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
\(\begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} l^{-1}_{11} & 0 & 0 & 0 \\ l^{-1}_{21} & l^{-1}_{22} & 0 & 0 \\ l^{-1}_{31} & l^{-1}_{32} & l^{-1}_{33} & 0 \\ l^{-1}_{41} & l^{-1}_{42} & l^{-1}_{43} & l^{-1}_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)
\(l_{11} l^{-1}_{11} = 1 \\ l_{21} l^{-1}_{11} + l_{22} l^{-1}_{21} = 0 \\ l_{31} l^{-1}_{11} + l_{32} l^{-1}_{21} + l_{33} l^{-1}_{31} = 0 \\ l_{41} l^{-1}_{11} + l_{42} l^{-1}_{21} + l_{43} l^{-1}_{31} + l_{44} l^{-1}_{41} = 0 \)
ie.\(l^{-1}_{11} = \frac{1}{l_{11}} \\ l^{-1}_{21} = -\frac{1}{l_{22}}(l_{21} l^{-1}_{11}) \\ l^{-1}_{31} = -\frac{1}{l_{33}}(l_{31} l^{-1}_{11} + l_{32} l^{-1}_{21}) \\ l^{-1}_{41} = -\frac{1}{l_{44}}(l_{41} l^{-1}_{11} + l_{42} l^{-1}_{21} + l_{43} l^{-1}_{31}) \)
\(l_{22} l^{-1}_{22} = 1 \\ l_{32} l^{-1}_{22} + l_{33} l^{-1}_{32} = 0 \\ l_{42} l^{-1}_{22} + l_{43} l^{-1}_{32} + l_{44} l^{-1}_{42} = 0 \)
ie.\(l^{-1}_{22} = \frac{1}{l_{22}} \\ l^{-1}_{32} = -\frac{1}{l_{33}}(l_{32} l^{-1}_{22}) \\ l^{-1}_{42} = -\frac{1}{l_{44}}(l_{42} l^{-1}_{22} + l_{43} l^{-1}_{32}) \)
\(l_{33} l^{-1}_{33} = 1 \\ l_{43} l^{-1}_{33} + l_{44} l^{-1}_{43} = 0 \)
ie.\(l^{-1}_{33} = \frac{1}{l_{33}} \\ l^{-1}_{43} = -\frac{1}{l_{44}}(l_{43} l^{-1}_{33}) \)
\(l_{44} l^{-1}_{44} = 1 \)
ie.\(l^{-1}_{44} = \frac{1}{l_{44}} \)
\(l^{-1}_{kk} = \frac{1}{l_{kk}} \\ l^{-1}_{ik} = -\frac{1}{l_{ii}} \sum_{j=k}^{i-1} l_{ij} l^{-1}_{jk} \\ (i=k+1,k+2,...,n) \)
\(\begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} u^{-1}_{11} & u^{-1}_{12} & u^{-1}_{13} & u^{-1}_{14} \\ 0 & u^{-1}_{22} & u^{-1}_{23} & u^{-1}_{24} \\ 0 & 0 & u^{-1}_{33} & u^{-1}_{34} \\ 0 & 0 & 0 & u^{-1}_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)
\(u^{-1}_{11} = 1 \)
\(u^{-1}_{12} + u_{12}u^{-1}_{22} = 0 \\ u^{-1}_{22} = 1 \)
ie.\(u^{-1}_{22} = 1 \\ u^{-1}_{12} = -(u_{12} u^{-1}_{22}) \)
\(u^{-1}_{13} + u_{12} u^{-1}_{23} + u_{13}u^{-1}_{33} = 0 \\ u^{-1}_{23} + u_{23} u^{-1}_{33} = 0 \\ u^{-1}_{33} = 1 \)
ie.\(u^{-1}_{33} = 1 \\ u^{-1}_{23} = -(u_{23} u^{-1}_{33}) \\ u^{-1}_{13} = -(u_{12} u^{-1}_{23} + u_{13}u^{-1}_{33}) \)
\(u^{-1}_{14} + u_{12} u^{-1}_{24} + u_{13} u^{-1}_{34} + u_{14}u^{-1}_{44} = 0 \\ u^{-1}_{24} + u_{23} u^{-1}_{34} + u_{24}u^{-1}_{44} = 0 \\ u^{-1}_{34} + u_{34}u^{-1}_{44} = 0 \\ u^{-1}_{44} = 1 \)
ie.\(u^{-1}_{44} = 1 \\ u^{-1}_{34} = -(u_{34} u^{-1}_{44}) \\ u^{-1}_{24} = -(u_{23} u^{-1}_{34} + u_{24} u^{-1}_{44}) \\ u^{-1}_{14} = -(u_{12} u^{-1}_{24} + u_{13} u^{-1}_{34} + u_{14} u^{-1}_{44}) \)
\(u^{-1}_{kk}=1 \\ u^{-1}_{ik} = -\sum_{j=i+1}^{k} u_{ij} u^{-1}_{jk} \\ (i=k-1,k-2,...,1) \)
/* * Copyright 2017 HONDOH Atsushi. */ package com.matarapi; import com.aparapi.Kernel; /** * Matrix calculate Kernel. * @author a.ho */ public abstract class MatKernel extends Kernel { /** * data array. */ protected float ary[]; /** * data offset. */ protected int[] offset; /** * matrix size. */ protected int[] matSize; /** * col size. */ protected int[] colSize; /** * Constructor. * @param data matrixes */ public MatKernel(IMatrix ...data) { int arySize = 0; for (IMatrix m : data) { arySize += m.getSize(); } ary = new float[arySize]; offset = new int[data.length]; matSize = new int[data.length]; colSize = new int[data.length]; int p = 0; int n = 0; for (IMatrix m : data) { offset[n] = p; matSize[n] = m.getRowSize() * m.getColSize(); colSize[n] = m.getColSize(); p = m.accept(ary, p); n += 1; } } /** * get Matrix. * @param no matrix number * @return Matrix */ public Matrix getMat(final int no) { return new Matrix(ary, offset[no], matSize[no] / colSize[no], colSize[no]); } /** * add. * <pre> * in1 in2 out * |a b c| |A B C| |a+A b+B c+C| * |d e f| + |D E F| -> |d+D e+E f+F| * |g h i| |G H I| |g+G h+H i+I| * * paralles size : matrix size (= row * col) * </pre> * @param in1 input1 * @param in2 input2 * @param out output */ protected void matAdd(int in1, int in2, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int p0 = offset[out] + i; int p1 = offset[in1] + i; int p2 = offset[in2] + i; ary[p0] = ary[p1] + ary[p2]; } /** * sub. * <pre> * in1 in2 out * |a b c| |A B C| |a-A b-B c-C| * |d e f| - |D E F| -> |d-D e-E f-F| * |g h i| |G H I| |g-G h-H i-I| * * paralles size : matrix size (= row * col) * </pre> * @param in1 input1 * @param in2 input2 * @param out output */ protected void matSub(int in1, int in2, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int p0 = offset[out] + i; int p1 = offset[in1] + i; int p2 = offset[in2] + i; ary[p0] = ary[p1] - ary[p2]; } /** * ReLU (Rectified linear unit). * <pre> * An simple activation function for neural network. * This is enough complicity for machine learning. * That's differentiate is so simple! * * ReLu(x) = max(0,x) * = 0 (x < 0) * x (x >= 0) * * cf. * ReLu-1(x) = 0 (x < 0) * 1 (x >= 0) * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output */ protected void matReLu(int in,int out) { int i = getGlobalId(); if (i > matSize[out]) return; int p0 = offset[out] + i; int p1 = offset[in] + i; ary[p0] = max(ary[p1],0.0f); } /** * differentiate of ReLU (Rectified linear unit). * <pre> * an simple activation function for neural network. * * ReLu-1(x) = 0 (x < 0) * 1 (x >= 0) * * cf. * ReLu(x) = max(0,x) * = 0 (x < 0) * x (x >= 0) * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output */ protected void matdReLu(int in,int out) { int i = getGlobalId(); if (i > matSize[out]) return; int p0 = offset[out] + i; int p1 = offset[in] + i; ary[p0] = ary[p1] < 0.0f ? 0.0f : 1.0f; } /** * multiple. * <pre> * in1 in2 out * |a b c| |A B C| |aA+bD+cG dA+eD+fG gA+hD+iG| * |d e f| X |D E F| -> |aB+bE+cH dB+eE+fH gB+hE+iH| * |g h i| |G H I| |aC+aF+aI dC+eF+fI gC+hF+iI| * * paralles size : matrix size (= row * col) * </pre> * @param in1 input1 * @param in2 input2 * @param out output */ protected void matMul(int in1, int in2, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int c0 = colSize[out]; int c1 = colSize[in1]; int c2 = colSize[in2]; int col = i % c0; int row = (i - col) / c0; int p0 = offset[out] + i; int p1 = offset[in1] + c1 * row; int p2 = offset[in2] + col; float mul = 0.0f; for (int cnt=0; cnt < c1; cnt++) { mul += ary[p1] * ary[p2]; p1 += 1; p2 += c2; } ary[p0] = mul; } /** * multiple LU. * <pre> * in1 in2 out * |a 0 0| |1 B C| |a aB aC | * |d e 0| X |0 1 F| -> |d dB+e dC+eF | * |g h i| |0 0 1| |g gB+h gC+hF+i | * * in1 and in2 is combined represnted in arg 'in' * |a B C| * |d e F| * |g h i| * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output */ protected void matMulLU(int in, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int size = colSize[in]; int col = i % size; int row = (i - col) / size; int p0 = offset[out] + i; int pL = offset[in] + size * row; int pU = offset[in] + col; float mul = 0.0f; int step = min(col, row); for (int cnt=0; cnt <= step; cnt++) { if (col == cnt) { mul += ary[pL]; } else { mul += ary[pL] * ary[pU]; } pL += 1; pU += size; } ary[p0] = mul; } /** * multiple UL. * <pre> * in1 in2 out * |1 B C| |a 0 0| |a+Bd+Cg Be+Ch Ci | * |0 1 F| X |d e 0| -> |d+Fg e+FH Fi | * |0 0 1| |g h i| |g h i | * * in1 and in2 is combined represnted in arg 'in' * |a B C| * |d e F| * |g h i| * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output */ protected void matMulUL(int in, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int size = colSize[in]; int col = i % size; int row = (i - col) / size; int skip = max(row, col); int p0 = offset[out] + i; int pU = offset[in] + size * row + skip; int pL = offset[in] + col + size * skip; float mul = 0.0f; for (int cnt=skip; cnt < size; cnt++) { if (row == cnt) { mul += ary[pL]; } else { mul += ary[pU] * ary[pL]; } pU += 1; pL += size; } ary[p0] = mul; } /** * ivert Upper triangle matrix. * <pre> * in out E * |1 b c| |1 B C| |1 0 0| * |0 1 f| X |0 1 F| -> |0 1 0| * |0 0 1| |0 0 1| |0 0 1| * * The Lower triangle of input matrix, includes diagonal components '1' * , is ignored. * * paralles size : col size * </pre> * @param in input * @param out output */ protected void matInvU(int in, int out) { int k = getGlobalId(); int size = colSize[in]; if (k > size) { return; } int p0 = offset[out] + size * k + k; ary[p0] = 1.0f; float inv; int pU, pUi; for (int i = k - 1; i >= 0; i--) { inv = 0.0f; int j = i + 1; pU = offset[in] + size * i + j; // U[i][j] pUi = offset[out] + size * j + k; // U-1[j][k] for (; j <= k; j++) { inv -= ary[pU] * ary[pUi]; pU += 1; // move right pUi += size; // move down } p0 -= size; // move up ary[p0] = inv; } } /** * ivert Lower triangle matrix. * <pre> * in out E * |a 0 0| |A 0 0| |1 0 0| * |d e 0| X |D E 0| -> |0 1 0| * |g h i| |G H I| |0 0 1| * * The Upper triangle of input matrix is ignored. * * paralles size : col size * </pre> * @param in input * @param out output */ protected void matInvL(int in, int out) { int k = getGlobalId(); int size = colSize[in]; if (k > size) { return; } int p0 = offset[out] + size * k + k; int p1 = offset[in] + size * k + k; ary[p0] = 1.0f / ary[p1]; float inv; int pL, pLi; for (int i = k + 1; i < size; i++) { inv = 0.0f; int j = k; pL = offset[in] + size * i + j; // L[i][j] pLi = offset[out] + size * j + j; // L-1[j][k] for (; j <= i-1; j++) { inv -= ary[pL] * ary[pLi]; pL += 1; // move right pLi += size; // move down } p0 += size; // move down // L-1[i][k] = - (¦²(L[i][j]*L-1[j][k]) / L[i][i] ary[p0] = inv / ary[offset[in] + size * i + i]; } } /** * calculate hadamard product. * <pre> * in1 in2 out * |a b c| |A B C| |aA bB cC| * |d e f| O |D E F| -> |dD eE fF| * |g h i| |G H I| |gG hH iI| * </pre> * @param in1 input1 * @param in2 input2 * @param out output */ protected void matHmul(int in1, int in2, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int p0 = offset[out] + i; int p1 = offset[in1] + i; int p2 = offset[in2] + i; ary[p0] = ary[p1] * ary[p2]; } /** * calculate transpose matrix. * <pre> * in out * |a b c| = |a d| * |d e f| |b e| * |c f| * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output */ protected void matTranspose(int in, int out) { int i = getGlobalId(); if (i > matSize[out]) return; int c0 = colSize[out]; int c1 = colSize[in]; int col = i % c1; int row = (i - col) / c1; int p1 = offset[in] + i; int p0 = offset[out] + row + col * c0; // col <-> row ary[p0] = ary[p1]; } /** * sort columns. * <pre> * in out * |a b c| = |b a c| * |d e f| |e d f| * |g h i| |h g i| * * order = [1,0,2] means: * Col 0 -> Col 1 * Col 1 -> Col 0 * Col 2 -> Col 2 * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output * @param order sort order */ protected void matSortCol(int in, int out, int[] order) { int i = getGlobalId(); if (i > matSize[out]) return; int c0 = colSize[out]; int c1 = colSize[in]; int col = i % c1; int row = (i - col) / c1; int p0 = offset[out] + order[col] + c0 * row; int p1 = offset[in] + i; ary[p0] = ary[p1]; } /** * sort rows. * <pre> * in out * |a b c| = |d e f| * |d e f| |a b c| * |g h i| |g h i| * * order = [1,0,2] means: * Row 0 -> Row 1 * Row 1 -> Row 0 * Row 2 -> Row 2 * * paralles size : matrix size (= row * col) * </pre> * @param in input * @param out output * @param order sort order */ protected void matSortRow(int in, int out, int[] order) { int i = getGlobalId(); if (i > matSize[out]) return; int c0 = colSize[out]; int c1 = colSize[in]; int col = i % c1; int row = (i - col) / c1; int p0 = offset[out] + col + c0 * order[row]; int p1 = offset[in] + i; ary[p0] = ary[p1]; } @Override public String toString() { StringBuilder sb = new StringBuilder(); sb.append("ary:\n"); for (float val : ary) { sb.append(String.format("%.2f ", val)); } sb.append("\noffset:\n"); for (int val : offset) { sb.append(String.format("%d ", val)); } sb.append("\nrow size:\n"); int cnt = 0; for (int val : matSize) { sb.append(String.format("%d ", val / colSize[cnt++])); } sb.append("\ncol size:\n"); for (int val : colSize) { sb.append(String.format("%d ", val)); } return sb.toString(); } }
/* * Copyright 2017 HONDOH Atsushi. */ package com.matarapi.test; import com.matarapi.IMatrix; import com.matarapi.MatKernel; import com.matarapi.Matrix; import com.matarapi.MockMatrix; import static com.matarapi.test.MatAssert.assertMatrix; import org.junit.Test; import static com.matarapi.test.MatAssert.assertMatrixInvert; /** * Simple calculation test. * * @author atsushi */ public class MatrixTest { @Test public void testLU() throws Exception { Matrix m0 = new Matrix(new float[][]{ {2.0f, 5.0f, 7.0f, 8.0f}, {4.0f, 13.0f, 20.0f, 25.0f}, {8.0f, 29.0f, 50.0f, 71.0f}, {10.0f, 34.0f, 78.0f, 98.0f}, }); // m0 will be over-writed to lu-matrix Matrix original = m0.copyAll(); // m0 => LU // oreder is row order for LU dividing final int[] order = m0.toLU(); // L U => mul IMatrix work = new MockMatrix(4, 4); MatKernel kernel = new MatKernel(m0, work) { @Override public void run() { int pass = getPassId(); if (0 == pass) { matMulLU(0, 1); } else if (1 == pass) { // revert pivoting matSortRow(1, 0, order); } } }; kernel.execute(m0.getSize(),2); Matrix mul = kernel.getMat(0); // System.out.println(mul.toString()); // LU expects m0 assertMatrix(original, mul, 1.0E-6f); } @Test public void testInvert() throws Exception { Matrix m0 = new Matrix(new float[][]{ {2.0f, 5.0f, 7.0f, 8.0f}, {4.0f, 13.0f, 20.0f, 25.0f}, {8.0f, 29.0f, 50.0f, 71.0f}, {10.0f, 34.0f, 78.0f, 98.0f}, }); // m0 will be over-writed to lu-matrix Matrix original = m0.copyAll(); // m0 => LU with pivoting // oreder is row order for LU dividing final int[] order = m0.toLU(); // L U => mul IMatrix work = new MockMatrix(4, 4); MatKernel kernel = new MatKernel(m0, work) { @Override public void run() { int pass = getPassId(); // Can't use switch statement on Aparapi. if (0 == pass) { // U => U(-1) matInvU(0,1); // L => L(-1) matInvL(0,1); } else if (1 == pass) { // U(-1) L(-1) => A(-1) matMulUL(1,0); } else if (2 == pass) { // revert pivoting matSortCol(0,1,order); } } }; kernel.execute(m0.getSize(), 3); Matrix invert = kernel.getMat(1); // System.out.println(expect); // System.out.println(invert); // precision is not well, but its ok. // It's enough for the machine learning. // We don't use GPGPU for Rocket Science. assertMatrixInvert(original, invert, 1.0E-5f); } }