1// Copyright ©2015 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4// 5// Some of the loop unrolling code is copied from: 6// http://golang.org/src/math/big/arith_amd64.s 7// which is distributed under these terms: 8// 9// Copyright (c) 2012 The Go Authors. All rights reserved. 10// 11// Redistribution and use in source and binary forms, with or without 12// modification, are permitted provided that the following conditions are 13// met: 14// 15// * Redistributions of source code must retain the above copyright 16// notice, this list of conditions and the following disclaimer. 17// * Redistributions in binary form must reproduce the above 18// copyright notice, this list of conditions and the following disclaimer 19// in the documentation and/or other materials provided with the 20// distribution. 21// * Neither the name of Google Inc. nor the names of its 22// contributors may be used to endorse or promote products derived from 23// this software without specific prior written permission. 24// 25// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 26// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 27// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 28// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 29// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 30// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 31// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 32// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 33// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 34// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 35// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 37//+build !noasm,!appengine,!safe 38 39#include "textflag.h" 40 41#define X_PTR SI 42#define Y_PTR DX 43#define DST_PTR DI 44#define IDX AX 45#define LEN CX 46#define TAIL BX 47#define ALPHA X0 48#define ALPHA_2 X1 49 50// func AxpyUnitaryTo(dst []float64, alpha float64, x, y []float64) 51TEXT ·AxpyUnitaryTo(SB), NOSPLIT, $0 52 MOVQ dst_base+0(FP), DST_PTR // DST_PTR := &dst 53 MOVQ x_base+32(FP), X_PTR // X_PTR := &x 54 MOVQ y_base+56(FP), Y_PTR // Y_PTR := &y 55 MOVQ x_len+40(FP), LEN // LEN = min( len(x), len(y), len(dst) ) 56 CMPQ y_len+64(FP), LEN 57 CMOVQLE y_len+64(FP), LEN 58 CMPQ dst_len+8(FP), LEN 59 CMOVQLE dst_len+8(FP), LEN 60 61 CMPQ LEN, $0 62 JE end // if LEN == 0 { return } 63 64 XORQ IDX, IDX // IDX = 0 65 MOVSD alpha+24(FP), ALPHA 66 SHUFPD $0, ALPHA, ALPHA // ALPHA := { alpha, alpha } 67 MOVQ Y_PTR, TAIL // Check memory alignment 68 ANDQ $15, TAIL // TAIL = &y % 16 69 JZ no_trim // if TAIL == 0 { goto no_trim } 70 71 // Align on 16-byte boundary 72 MOVSD (X_PTR), X2 // X2 := x[0] 73 MULSD ALPHA, X2 // X2 *= a 74 ADDSD (Y_PTR), X2 // X2 += y[0] 75 MOVSD X2, (DST_PTR) // y[0] = X2 76 INCQ IDX // i++ 77 DECQ LEN // LEN-- 78 JZ end // if LEN == 0 { return } 79 80no_trim: 81 MOVQ LEN, TAIL 82 ANDQ $7, TAIL // TAIL := n % 8 83 SHRQ $3, LEN // LEN = floor( n / 8 ) 84 JZ tail_start // if LEN == 0 { goto tail_start } 85 86 MOVUPS ALPHA, ALPHA_2 // ALPHA_2 := ALPHA for pipelining 87 88loop: // do { 89 // y[i] += alpha * x[i] unrolled 8x. 90 MOVUPS (X_PTR)(IDX*8), X2 // X_i = x[i] 91 MOVUPS 16(X_PTR)(IDX*8), X3 92 MOVUPS 32(X_PTR)(IDX*8), X4 93 MOVUPS 48(X_PTR)(IDX*8), X5 94 95 MULPD ALPHA, X2 // X_i *= alpha 96 MULPD ALPHA_2, X3 97 MULPD ALPHA, X4 98 MULPD ALPHA_2, X5 99 100 ADDPD (Y_PTR)(IDX*8), X2 // X_i += y[i] 101 ADDPD 16(Y_PTR)(IDX*8), X3 102 ADDPD 32(Y_PTR)(IDX*8), X4 103 ADDPD 48(Y_PTR)(IDX*8), X5 104 105 MOVUPS X2, (DST_PTR)(IDX*8) // y[i] = X_i 106 MOVUPS X3, 16(DST_PTR)(IDX*8) 107 MOVUPS X4, 32(DST_PTR)(IDX*8) 108 MOVUPS X5, 48(DST_PTR)(IDX*8) 109 110 ADDQ $8, IDX // i += 8 111 DECQ LEN 112 JNZ loop // } while --LEN > 0 113 CMPQ TAIL, $0 // if TAIL == 0 { return } 114 JE end 115 116tail_start: // Reset loop registers 117 MOVQ TAIL, LEN // Loop counter: LEN = TAIL 118 SHRQ $1, LEN // LEN = floor( TAIL / 2 ) 119 JZ tail_one // if LEN == 0 { goto tail } 120 121tail_two: // do { 122 MOVUPS (X_PTR)(IDX*8), X2 // X2 = x[i] 123 MULPD ALPHA, X2 // X2 *= alpha 124 ADDPD (Y_PTR)(IDX*8), X2 // X2 += y[i] 125 MOVUPS X2, (DST_PTR)(IDX*8) // y[i] = X2 126 ADDQ $2, IDX // i += 2 127 DECQ LEN 128 JNZ tail_two // } while --LEN > 0 129 130 ANDQ $1, TAIL 131 JZ end // if TAIL == 0 { goto end } 132 133tail_one: 134 MOVSD (X_PTR)(IDX*8), X2 // X2 = x[i] 135 MULSD ALPHA, X2 // X2 *= a 136 ADDSD (Y_PTR)(IDX*8), X2 // X2 += y[i] 137 MOVSD X2, (DST_PTR)(IDX*8) // y[i] = X2 138 139end: 140 RET 141