subroutine Dsolve_Mrhs(A,b,trans)
real(8),dimension(:,:),intent(in) :: A
real(8),dimension(:,:),intent(inout) :: b
character(len=1),optional :: trans
character(len=1) :: trans_
integer :: m,n,nrhs,lda,ldb
integer :: info
integer,dimension(:),allocatable :: ipvt
trans_="N";if(present(trans))trans_=trans
lda = max(1,size(A,1))
ldb = max(1,size(B,1))
m = size(A,1)
n = size(A,2)
nrhs= size(B,2)
allocate(ipvt(min(m,n)))
call dgetrf(m,n,A,lda,ipvt,info)
if(info/=0)stop "Error MATRIX/d_mat_solve_linear_system: dgetrf"
call dgetrs(trans_,n,nrhs,A,lda,ipvt,b,ldb,info)
if(info/=0)stop "Error MATRIX/d_mat_solve_linear_system: dgetrs"
deallocate(ipvt)
end subroutine dsolve_Mrhs