1/* fflas/fflas_ftrsv.inl 2 * Copyright (C) 2005 Clement Pernet 3 * 4 * Written by Clement Pernet <Clement.Pernet@imag.fr> 5 * 6 * 7 * ========LICENCE======== 8 * This file is part of the library FFLAS-FFPACK. 9 * 10 * FFLAS-FFPACK is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 23 * ========LICENCE======== 24 *. 25 */ 26#ifndef __FFLASFFPACK_ftrsv_INL 27#define __FFLASFFPACK_ftrsv_INL 28 29namespace FFLAS { 30 //--------------------------------------------------------------------- 31 // ftrsv: TRiangular System solve with vector 32 // Computes X <- op(A^-1).X 33 // size of X is m 34 //--------------------------------------------------------------------- 35 template<class Field> 36 inline void 37 ftrsv (const Field& F, const FFLAS_UPLO Uplo, 38 const FFLAS_TRANSPOSE TransA, const FFLAS_DIAG Diag, 39 const size_t N,typename Field::ConstElement_ptr A, size_t lda, 40 typename Field::Element_ptr X, int incX) 41 { 42 43 typename Field::Element_ptr Xi; 44 typename Field::ConstElement_ptr Ai; 45 if ( Uplo == FflasLower ){ 46 if ( TransA == FflasTrans){ 47 Ai = A+(N-1)*(lda+1); // bottom right entry of A 48 Xi = X+(int)(N-1)*incX; 49 for(size_t i=0; i<N; Ai-=(lda+1), Xi-=incX, ++i) { 50 if (i) 51 F.subin(*Xi, fdot(F, i, Ai+lda, lda, Xi+incX, incX)); 52 if ( Diag==FflasNonUnit ) 53 F.divin(*Xi,*Ai); 54 } 55 } // End FflasTrans 56 else{ 57 Ai = A; 58 Xi = X; 59 for(size_t i=0 ; i<N; Ai+=lda,Xi+=incX, ++i ){ 60 F.subin (*Xi, fdot (F, i, Ai, 1, X, incX)); 61 if ( Diag==FflasNonUnit ) 62 F.divin(*Xi,*(Ai+i)); 63 } 64 } 65 } // End EFflasLower 66 else{ 67 if ( TransA == FflasTrans){ 68 Ai = A; 69 Xi = X; 70 for(size_t i=0; i<N; ++Ai,Xi+=incX, ++i) { 71 F.subin(*Xi, fdot(F, i, Ai, lda, X, incX)); 72 if ( Diag == FflasNonUnit ) 73 F.divin(*Xi, *(Ai+i*lda)); 74 } 75 76 } // End FflasTrans 77 else{ 78 Ai = A+(lda+1)*(N-1); 79 Xi = X+incX*(int)(N-1); 80 for(size_t i=0; Xi>=X; Ai-=lda+1,Xi-=incX, ++i ){ 81 if (i) 82 F.subin (*Xi, fdot (F, i, Ai+1, 1, Xi+incX, incX)); 83 if ( Diag==FflasNonUnit ) 84 F.divin(*Xi,*Ai); 85 } 86 } 87 } 88 } 89 90} 91#endif // __FFLASFFPACK_ftrsv_INL 92/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 93// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 94